Life-long telomere attrition predicts health and lifespan in a large mammal
Luise A. Seeker
Thursday, April 20, 2017
PART 1: THE MANUSCRIPT
Short running title
Telomere length attrition predicts health span
Abstract
Telomere length measured in blood cells is predictive of subsequent adult health and survival across a range of vertebrate species. However, we currently do not know whether such associations result from among-individual differences in telomere length determined genetically or by environmental factors early in life, or from differences in the rate of telomere attrition over the course of life. Here, we measured relative leukocyte telomere length (RLTL) multiple times across the entire lifespan of dairy cattle in a research population that is closely monitored for health and milk production and where individuals are only culled in response to health issues and less due to poor milk production than on purely commercial farms. Our results clearly show that the average amount of telomere attrition over an individual’s life, not their average or early life telomere length predicted when an individual was culled. Within-individual telomere length attrition could reflect environmental or physiological insults which may accumulate to predict individual health-span. We also show that animals with more telomere attrition in their first year of life were culled at a younger age, indicating that early life stressors may have a prolonged effect on adult life.
Keywords:
Telomeres, telomere length, longitudinal, lifespan, health-span, qPCR
Introduction
Telomeres are repetitive DNA sequences that cap the ends of eukaryote linear chromosomes (Blackburn & Gall, 1978; De Lange, 2005). They shorten with the number of cell divisions in vitro as well as in response to oxidative stress and critically short telomeres trigger a DNA damage response that leads to replicative senescence or apoptosis (Harley, Futcher, & Greider, 1990; Olovnikov, 1973; Watson, 1972). In the last decade or so, measures of average telomere length (TL) taken from blood samples have emerged as an exciting biomarker of health across disciplines including biomedicine, epidemiology, ecology and evolutionary biology (Aviv & Shay, 2018; Harrington & Pucci, 2018; Wilbourn et al., 2018). Considerable among- and within-individual variation in TL has been observed, with a general pattern of rapid telomere attrition during early life and a plateau or slower decline thereafter (Aubert & Lansdorp, 2008; Baerlocher, Rice, Vulto, & Lansdorp, 2007). Both genetic and environmental factors, particularly those associated with physiological stress, predict TL in humans and other vertebrates (Angelier, Costantini, Blévin, & Chastel, 2017; Muhammad Asghar et al., 2015; Cherkas et al., 2006; Dugdale, Richardson, & Richardson, 2018; Epel et al., 2004). TL has also been repeatedly associated with health outcomes and subsequent survival in a variety of species, particularly humans and birds (Boonekamp, Simons, Hemerik, & Verhulst, 2013; Wilbourn et al., 2018) and also experimentally elongated TL in mice were associated with a survival advantage (Muñoz-lorente, Cano-martin, & Blasco, 2019). However, a major outstanding question remains to what degree associations between TL and health arise from constitutive differences in TL among individuals set by genes or early life conditions, or from the pattern of within-individual change in TL across individuals’ lives. Estimates of the individual consistency of TL over time in both human and avian literature vary considerably among studies. Some studies report very high correlations within individuals across follow-up measurements (r > 0.8; (Benetos et al., 2013; Boonekamp, Mulder, Salomons, Dijkstra, & Verhulst, 2014)) and high heritability of TL (>0.7; (Atema et al., 2015; Broer et al., 2013)), and have shown in humans that the rank order in TL among individuals remains relatively consistent during adult life (Benetos et al., 2013). This implies most of the variation in blood cell TL occurs at the among-individual level and is predominantly determined by genetics and early-life environment (Atema et al., 2015; Benetos et al., 2013). In stark contrast to this, a growing body of literature studying both humans and non-human vertebrates reports much lower individual repeatability and heritability (Dugdale & Richardson, 2018; Voillemot et al., 2012). These studies demonstrate very high levels of within-individual variation in TL and show that changes in TL across consecutive measurements are highly dynamic (Fairlie et al., 2015; Farzaneh-Far et al., 2010). Although on average TL attrition over time tends to be the norm, a growing number of studies have found that a substantial proportion of individuals show telomere elongation over time and it is now clear this cannot be ascribed solely to measurement error (Bateson & Nettle, 2016; Spurgin et al., 2018). A recent cross-sectional study has shown that amongst different avian and mammal species, those with a fast telomere attrition rate have a shorter life span but without investigating if telomere attrition rate within individuals is predictive of their life span within the population as well (Whittemore, Vera, Martínez-nevado, Sanpera, & Blasco, 2019). To date, few studies have directly compared the power of an individual’s average or early life TL versus the pattern of within-individual change in TL to explain variation in measured associated with health and fitness. In Seychelle warblers and Alpine swifts both telomere length and telomere attrition rate were associated with survival (Barrett, Burke, Hammers, Komdeur, & Richardson, 2013; Bize, Criscuolo, Metcalfe, Nasir, & Monaghan, 2009) while in jackdaw nestlings telomere length was not predictive of lifespan but early life telomere attrition- accelerated by early life stress- was (Boonekamp, Mulder, Salomons, Dijkstra, & Verhulst, 2014). Contrary, in a study on elderly humans, there was an association between telomere length and longevity, but not between telomere attrition and life span (Yuan, Kronström, Hellenius, Cederholm, & Xu, 2018). To our knowledge no studies so far has investigated the relationship between telomere attrition and life span with data spanning the entire typical lifespan of the studied individuals. Here, we use an exceptionally detailed longitudinal study of dairy cattle to show that a life-long tendency for greater telomere attrition, rather than average TL, predicts productive lifespan, a metric closely linked to long-term health.
We used samples and data collected as part of the long-term study of Holstein Friesian dairy cattle kept at the Crichton Royal Research Farm in Dumfries, Scotland (Veerkamp, Simm, & Oldham, 1994). This herd, consisting of around 200 milking cows plus their calves and replacement heifers, has been regularly monitored since 1973 for a broad range of measurements, such as body weight, feed intake, signs of disease (health events), milk yield, productive lifespan and reasons for culling (Veerkamp et al., 1994). The cows are separated into two distinct genetic groups (a group that has been strongly selected for high milk fat and protein yield, and a control group) and are randomly allocated to two different diets (a high forage, low energy diet vs. a low forage, high energy diet group) (Veerkamp et al., 1994). Blood samples are routinely collected from members of the herd, initially within 15 days after birth and then approximately annually thereafter. Productive lifespan, which is the age of an individual at culling, is recorded for every animal together with a reason for culling. In this and other commercial livestock herds, culling is a response to diseases such as reproductive problems, lameness and mastitis. Although productive lifespan in livestock is clearly not a comparable measure to longevity in study of humans or wild vertebrates, we argue that because age at culling is tightly linked to health, it can be considered as a proxy for health-span, over and above its obvious importance as a metric of agricultural productivity. We have previously reported for this herd that TL is moderately repeatable and heritable, but that individual slopes relating TL to age vary considerably among individuals (Seeker et al., 2018; Seeker et al., 2018). We have also reported that TL at the age of one year predicts productive lifespan in this population of dairy cattle (Seeker et al., 2018). Here, we test the relationships between lifetime average TL as well as lifetime change in TL with productive lifespan.
Materials and Methods
Animal population and data collection
At the Crichton Royal Farm, 200 milking cows plus their calves and replacement heifers are kept at any time. One half of the milking cows belong to a genetic line that has been selected for high milk protein and fat yield (S), while the other half is deliberately maintained on a UK average productivity level (C). Calves and heifers of both genetic lines are kept together. After first calving all cows are randomly allocated to a high forage (HF) or low forage (LF) diet. The LF diet is energy richer than the HF diet and whilst the LF cows are housed continuously, the HF cows graze over the summer months. All cows are milked three times daily and milk yield is recorded. In the present study, these measurements were used to calculate an average milk production in kg per cow including all started lactations and it is referred to this below as “average lifetime milk production” (Figure S1). Every day cows leave the milking parlour over a pressure plate which detects signs of lameness. Behaviour and health events are documented after visual detection by farm workers (Figure S2). At the end of the animal’s life its productive lifespan and a reason for culling are recorded. Productive lifespan is the time from birth to culling in days and is a proxy for the health span of the animal, because all animals that remain healthy enough to generate profit for the farmer remain in the herd. The most frequent reasons for culling were reproductive problems, mastitis, lameness and injuries caused by accidents (Figure S3). Routine blood sampling takes place initially shortly after birth and then annually in spring. Because of this sampling routine and because calves are being born all year round, age at sampling and sampling intervals vary for adult animals (Figure S1). Further information about the animal population can be found in Supplementary File 2.
DNA extraction and qPCR
DNA from whole blood samples was extracted with the DNeasy Blood and Tissue spin column kit (QIAGEN) and telomere length was measured by qPCR as previously described (Seeker et al., 2016; Seeker et al., 2018; Seeker et al., 2018; Seeker, Holland, Psifidi, Banos, & Nussey, 2015). The repeatability of the assay (see Supplementary File 2 for how repeatability was calculated) was 80% and therefore delivers robust results for further interpretation (Nettle, Seeker, Nussey, Froy, & Bateson, 2019). A full description of our DNA extraction and qPCR protocols including quality control steps can also be found in Supplementary File 2.
Statistical analysis
We have shown before that our RLTL data are significantly affected by qPCR plate and qPCR row (Seeker et al., 2016; Seeker et al., 2018). To account for those known sources of measurement error, we used the residuals of a linear model that corrected all RLTL measurements for qPCR plate and row, by fitting plate and row as fixed factors in the model. These residual RLTL measures were used in all subsequent calculations and models of telomere dynamics. RLTL change was calculated as the difference between two subsequent adjusted RLTL measurements within individuals (RLTL change = RLTLt - RLTLt-1). To investigate which factors affect the direction and amount of RLTL change, a linear mixed model was fitted with RLTL change as response variable and animal identity as random effect. To capture the experimental design, feed group and genetic group were added as fixed effects to all models. Also, the time difference between consecutive samplings in days was fitted as a covariate, becasue sampling intervals varied within and between animals. The following additional factors were tested for their relationship with RLTL change: birth year of the animal (fixed effect), year of sampling (fixed effect), age at sampling (at time t), and the occurrence of a health event within two weeks before or after sampling (at time t).
Becasue of our study design, the first two sample years only included RLTL change in young animals which has previously been shown to follow a different pattern. To exclude that a year effect was driven by the difference in juvenile and adult RLTL dynamics, the full model was repeated while excluding the first two sample years.
Birth year and the occurrence of a heath event were not statistically significant (p>0.05) and therefore eliminated from the model. In the reduced model, age was first fitted as a covariate (age in years) than as a fixed effect (factor with two levels: younger vs. older than 1 year) and model fit was compared using AIC values. Because we wanted to investigate if milk productivity was associated with RLTL change, the best-fitting model of RLTL change was repeated for 918 change measurements of 253 animals for which milk productivity measurements were available and average lifetime milk production re-scaled by dividing it by 1000 was fitted as additional covariate.
To investigate the association between RLTL change and productive lifespan we first focussed on RLTL change within the first year of life by only considering two RLTL measurements per animal: The first was taken shortly after birth and the second at the approximate age of one year. The change between those measurements was noramlly distributed and tested as an explanatory variable in a Cox proportional hazard model of productive lifespan.
Our data shows similar to other longitudinal studies with more than two measurements no constant attrition, maintenance or elongation, but more complex dynamics with short term changes in both directions. Example individual RLTL dynamics are shown for all cows with at least seven RLTL measurements in our dataset in figure S4. The complexity of telomere dynamics (which cannot be explained by measurement error alone) means that it is impossible to compare lifelong telomere dynamics by simply comparing slopes. To still be able to test how life-long telomere dynamics are associated with productive lifespan, we calculated three average metrics of RLTL dynamics over every individual’s lifetime: Firstly, the average of all their RLTL measures (“mean RLTL”) which investigates if animals with on average longer RLTL have a survival advantage, secondly, the average of all their RLTL change measures (“mean RLTL change”) and thirdly, the average of all their absolute RLTL change measures (“mean absolute RLTL change”). While mean RLTL change captures the direction and magnitude of RLTL changes, mean absolute RLTL change describes just the magnitude of change without considering its direction because we were interested to investigate if more change in either direction may be correlated with adverse effects. Figure S5 visualises the reasoning behind calculating these three measures of lifetime telomere length dynamics which are not surprisingly moderately correlated with one another (r ranged from -0.53 to 0.30, Figure S6). We used Cox proportional hazard models of productive lifespan that also included mean milk production as covariate and included mean RLTL, mean RLTL change and mean absolute RLTL change as explanatory variables first separately and then together to test their association with productive lifespan first separately, then while accounting for the effect of the other two measures. We have reported before that on average RLTL shortens in this population of dairy cattle within the first year of life while remaining relatively stable thereafter at the population level (Seeker et al., 2018; Seeker, et al., 2018). We were concerned that the initial RLTL attrition may contribute proportionally more to overall RLTL shortening in short-lived animals with less follow up samples which may explain a correlation of more telomere attrition with a shorter productive lifespan. To ensure the relationship was not entirely driven by this correlation, we tested all models again while excluding all measurements that were taken shortly after birth and then also while excluding all animals with less than three RLTL measurements. Lastly, we tested if the previously reported effect of RLTL at the age of one year on productive lifespan (Seeker, Ilska, Psifidi, Wilbourn, Underwood, Fairlie, Holland, Froy, Salvo-Chirnside, et al., 2018) remained statistically significant when tested in a Cox proportional hazard model together with milk productivity and mean RLTL change. All statistical analyses were performed in R studio with R 3.1.3. (R Core Team, 2014). Mixed-effects models were implemented using the ‘lme4’ library (Bates, Mächler, Bolker, & Walker, 2015), while Cox proportional hazard models were implemented using the library survival (Therneau, 2015) and figures were generated with the library ‘ggplot2’ (Wickham, 2009).
Results
We used blood samples collected between 2008 and 2014 to measure relative leukocyte telomere length (RLTL) by monoplex qPCR in 1,325 samples from 305 individuals. On average 4.3 (range: 2-8) telomere measurements were made of each individual, including the first measurement within 15 days of birth and a variable number of subsequent measures (Figure S7 A). RLTL measures, adjusted for qPCR plate and row (see Methods), were approximately normally distributed (Figure S7 B).
From these longitudinal data we calculated 1,020 measures of change in adjusted RLTL across subsequent measurements of individual animals. RLTL change also showed approximately a normal distribution with a mean statistically significantly smaller than zero (Figure S7 C; t = -5.844, df = 1019, p-value <0.001) indicating that telomere shortening overall predominated.
Animals varied in the amount and direction of RLTL change across consecutive measurements, with a fairly even proportion of individuals increasing (43.5%) and decreasing (56.5%) in RLTL over time.
Modelling
Factors affecting telomere change
We next ran a series of models to test whether known individual, genetic and environmental variables could explain variation in RLTL change. Only sample year was significant in these models; age in years, genetic group, feed group, birth year, the time difference between samplings in days, and the occurrence of a health event within two weeks of sampling were not significant (Table S1). However, when age was fitted a factor with two levels for animals younger or oler than 1 year, it became clear that younger anumals show significantly more telomere length attrition (Tables S2 and S3). Fitting age as a covariate (age in years) and not a two-level factor lead to a slightly better although not statistically significant model fit and was therefore realised in all following models (AICs: -601.82 vs. -600.99).
We have previously shown that, on average across this population, RLTL declined over the first year of life but showed no systematic change with age thereafter (Seeker, Ilska, Psifidi, Wilbourn, Underwood, Fairlie, Holland, Froy, Salvo-Chirnside, et al., 2018; Seeker, Ilska, Psifidi, Wilbourn, Underwood, Fairlie, Holland, Froy, Bagnall, et al., 2018). Consistent with this, we found that average RLTL change across consecutive measurements was only significantly negative (indicating a tendency for attrition over time) when the first measurement was made close to birth and the follow up measurement at the age of around 1 year (Figure 1A, estimate= -0.115, SE= 0.01, t= -11.512, p<0.001; Table S3).
Average lifetime milk productivity was not statistically associated with RLTL change, when tested in a subset of animals (253 animals with 918 RLTL change measurements) that had milk productivity measurements available (Table S2).
Consecutive RLTL measurements made on the same individual were overall moderately positively and significantly correlated (r=0.38, p<0.001; Figure 1B), supporting our previously reported moderate and significant individual repeatability of RLTL (Seeker, Ilska, Psifidi, Wilbourn, Underwood, Fairlie, Holland, Froy, Bagnall, et al., 2018). We next tested whether measured RLTL change predicted the productive lifespan of cattle. Of all dead cows (N=244) the vast majority (N=241) had survived to their first lactation, but there was considerable variation in productive lifespan beyond this point (Figure S7 D). We used a Cox proportional hazard model to test whether early life RLTL change between two samples, one taken shortly after birth and the next at an approximate age of one year, predicted productive lifespan and found that greater early life RLTL attrition was associated with a shorter productive lifespan (N= 291, coefficient= -1.141, SE= 0.391, p= 0.004, Figure 2, Table S4). When each measure of lifetime RLTL dynamics was tested in a separate cox proportional hazard model of productive lifespan, mean RLTL did not significantly predict productive lifespan (coefficient= 0.341, SE= 0.591, p=0.564, Figure 2B , Table S4), but both mean RLTL change (coefficient= -5.209, SE= 0.845, p<0.001; Figure 2D, Table S4) and mean absolute RLTL change (coefficient= 2.939, SE= 0.970, p=0.002; Figure 2C, Table S4) were significantly associated with productive lifespan. When all three measures of lifetime RLTL dynamics were included in the same model only mean RLTL change remained significant (coefficient = -4.758, SE= 1.018, p<0.001; Table S4). This implies that the relationship between productive lifespan and mean absolute RLTL change was largely due to covariance with mean RLTL change. Thus, individuals that experienced greater telomere attrition over their lifetimes had a shorter productive lifespan and direction of RLTL change (rather than simply absolute magnitude) was an important aspect of this relationship. To visualise the relationship between RLTL change measurements and productive lifespan using Kaplan-Meier plots (Figure 2), continuous RLTL measures were transformed to a discrete scale by grouping them into tertiles (Figure S8). Cox proportional hazard models based on these tertile groupings of RLTL measures showed similar results to those reported above (Figure 2). The relationship between mean RLTL change and productive lifespan was robust to the inclusion of milk production (a physiological stressor, which is positively associated with productive lifespan as cows with a low milk yield are more likely to be culled at a younger age) in the model (Table S5). If RLTL declines most within the first year of life (as Figure 1A indicates), the association between change in RLTL and productive lifespan may be driven by the fact that this initial decline contributes relatively more to estimates of mean RLTL change in shorter-lived individuals than in animals that have more follow up samples with more moderate change measures available. We therefore repeated the analysis excluding the early life RLTL measurements and found that mean RLTL change still predicted productive lifespan (N=253, coefficient= -5.056, SE= 1.315, p<0.001, Table S6). The relationship between mean RLTL change and productive lifespan also remained statistically significant, when animals with less than 3 samples were excluded from the analysis (N= 213, coefficient= -3.47, SE= 1.34, Wald test = 6.7 on 1 df, p= 0.01) and when RLTL at 1 year of age (which we previously found to predict productive lifespan (Seeker, Ilska, Psifidi, Wilbourn, Underwood, Fairlie, Holland, Froy, Salvo-Chirnside, et al., 2018)) was included in the models (Table S7).
Mean rate of RLTL change and mean rate of absolute RLTL change correlate well with mean RLTL change and mean absolute RLTL change. THerefore, they do not seem to add much to the analysis. Only mean RLTL, mean RLTL change and mean absolute RLTL change are used.
Discussion
We have presented, to our knowledge, the first demonstration in any vertebrate that lifelong variation in telomere attrition rather than variation in constitutive individual differences in average telomere length predict health outcomes and lifespan. While there is mounting evidence that telomere length predicts mortality, health and life history in humans as well as birds and non-human mammals (Bakaysa et al., 2007; Bichet et al., 2019; Boonekamp et al., 2013; Cawthon, Smith, O’Brien, Sivatchenko, & Kerber, 2003; Fairlie et al., 2015; Mark Haussmann, Winkler, & Vleck, 2005; Heidinger et al., 2012; Nettle, Monaghan, Boner, Gillespie, & Bateson, 2013; Voillemot et al., 2012; Wilbourn et al., 2018), very few studies have been able to accumulate long-term longitudinal data capable of differentiating the role of among- and within-individual variation in TL to such relationships. Our data support the contention that within-individual directional change over time in TL is more important than stable among-individual differences in predicting overall health. Individuals in our study which showed a greater propensity to shorten RLTL across sampling points were culled earlier, mainly due to health issues described above. There was no relationship between productive lifespan and an individual’s average RLTL. Future studies will show how well our results generalise to other systems as telomere biology is variable amongst species. Cattle telomere biology seems to be similar to other ruminants, horses, zebras, tapirs some whales and primates including humans in that they have relatively short telomeres and a tight regulation of telomerase expression (Gomes et al., 2011). If our results generalise to some of those other systems and contexts, they have important implications for the utility of TL as a biomarker of health and fitness, lending support to the idea that TL change is an indirect marker reflecting past physiological insults and stress rather than an indicator of constitutive or genetically-based robustness to life’s challenges. Our data also highlight the importance of collecting longitudinal telomere measurements, by showing that in some species it is within-individual change over time in TL that carries the important biological signal.
Our study animals varied considerably in the magnitude and direction of RLTL change across consecutive sampling points. This is in accordance with other longitudinal studies that have reported a wide variation in TL change, alongside evidence that a large proportion of individuals show telomere lengthening over time (Bize et al., 2009; Chen et al., 2011; Farzaneh-Far et al., 2010; Gardner et al., 2005; Huzen et al., 2014; Kark, Goldberger, Kimura, Sinnreich, & Aviv, 2012; Nordfjäll et al., 2009; Shalev et al., 2012; Steenstrup et al., 2013; Svenson et al., 2011). Whilst the vast majority of studies linking telomere length with health and fitness remain cross-sectional, a growing number of longitudinal studies have found that individuals that lose more TL across samples have reduced survival prospects (Barrett et al., 2013; Bize et al., 2009; Boonekamp et al., 2014; Goglin et al., 2016). However, these studies tend to involve two samples taken per individual or a follow-up period that represents only a small fraction of the total lifespan of the study species. A lifelong study of TL and lifespan in zebra finches found that TL measured in early life was a better predictor of subsequent survival than measures from later life, but the role of change in telomere length across sampling points was not investigated (Heidinger et al., 2012). Here we show that, despite TL being moderately consistent across the lifetimes of individuals, there remains considerable within-individual variation and the pattern of change in TL over an individual’s life is highly dynamic. Short-term environmental fluctuations impacting TL dynamics could be responsible and may impact individuals in different ways. In humans many lifestyle factors such as socio-economic status, psychological stress, smoking and body mass index are known to influence telomere length (Cherkas et al., 2008; Epel, Daubenmier, Moskowitz, Folkman, & Blackburn, 2009; Huzen et al., 2014; Kim et al., 2009; Shalev et al., 2012). In animal studies chronic infection and exposure to stressful situations have been linked to telomere shortening (Asghar, Hasselquist, Zehtindjiev, Westerdahl, & Bensch, 2015; Nettle et al., 2013; Reichert et al., 2014). Whilst we were unable to identify variables other than age which predicted within-individual TL change in our study system, individuals with a greater propensity to lose TL over time had shorter productive lives, implying changes in TL reflect important environmental or physiological variation linked to health. From a telomere biology perspective, it would be interesting to investigate changes not only in mean TL but also in its variance over time, preferentially using a technique like STELA or TESLA that allows the detection of critically short telomeres as well (Baird, Rowson, Wynford-Thomas, & Kipling, 2003; Lai et al., 2017). However, if the focus is to find good markers for animal breeding and health those techniques are not economically practical.
Our results have important implications for the application of TL as a biomarker within animal breeding and health. Productive lifespan describes the time from birth to culling and differs from measures of “longevity” or “lifespan” in human or wild animal studies because the animals are culled for reasons relating mainly to health and productivity. From an animal production perspective, a biomarker which could predict the productive lifespan of livestock from a relatively early age would represent an exciting target phenotype for use in selective breeding if heritable or for use in animal welfare monitoring if influenced by the environment. In the present study we showed that RLTL attrition predicted productive lifespan. The fact that both early life and lifelong RLTL shortening predicts earlier culling does support recent suggestions that TL may be useful as a biomarker of stress and welfare in livestock (Melissa Bateson, 2016; Brown, Dechow, Liu, Harvatine, & Ott, 2012). Many previous studies in birds and humans have found that experience of past stress is associated with shorter telomeres (Angelier, Vleck, Holberton, & Marra, 2013; Asghar et al., 2015; Cherkas et al., 2006; Entringer et al., 2011; Entringer et al., 2013; Epel et al., 2004; Haussmann, Longenecker, Marchetto, Juliano, & Bowden, 2012; Nettle et al., 2017). In our study, cattle showing RLTL attrition may be those that have experienced some form of recent immediate stress or minor infection. Repeated bouts of such stress, reflected by higher average RLTL attrition through life, may result in the development of health issues that lead to culling. Whilst the results of the present study suggest that longitudinal RLTL change itself may be a useful biomarker of health and welfare in livestock, a key area for further study is the causes of variation in RLTL change, as these may help us to understand sources of physiological stress for domestic animals and help improve animal welfare and productivity.
Acknowledgements
We thank everyone involved in collecting blood samples on the Crichton farm, particularly Hannah Dykes and Isla McCubbin. We also thank the SynthSys facility at the University of Edinburgh, particularly Eliane Salvo-Chirnside, for their support with the liquid handling robot. Our work was funded by the Biotechnology and Biological Sciences Research Council (BB/L007312/1; http://www.bbsrc.ac.uk/) awarded to GB. LS received a studentship from Scotland’s Rural College (http://www.sruc.ac.uk/).
References
Angelier, F., Costantini, D., Blévin, P., & Chastel, O. (2017). Do glucocorticoids mediate the link between environmental conditions and telomere dynamics in wild vertebrates? A review. General and Comparative Endocrinology, 256, 99–111. https://doi.org/10.1016/j.ygcen.2017.07.007
Angelier, F., Vleck, C. M., Holberton, R. L., & Marra, P. P. (2013). Telomere length, non-breeding habitat and return rate in male American redstarts. Functional Ecology, 27(2), 342–350. https://doi.org/10.1111/1365-2435.12041
Asghar, M, Hasselquist, D., Zehtindjiev, P., Westerdahl, H., & Bensch, S. (2015). Hidden costs of infection: Chronic malaria accelerates telomere degradation and senescence in wild birds, 347(6220), 9–12.
Asghar, Muhammad, Bensch, S., Tarka, M., Hansson, B., Hasselquist, D., & Asghar, M. (2015). Maternal and genetic factors determine early life telomere length. Proceedings of the Royal Society B.
Atema, E., Mulder, E., Dugdale, H. L., Briga, M., van Noordwijk, A. J., & Verhulst, S. (2015). Heritability of telomere length in the Zebra Finch. Journal of Ornithology, 156(4), 1113–1123. https://doi.org/10.1007/s10336-015-1212-7
Aubert, G., & Lansdorp, P. M. (2008). Telomeres and aging. Physiological Reviews, 88(2), 557–579. https://doi.org/10.1152/physrev.00026.2007
Aviv, A., & Shay, J. W. (2018). Reflections on telomere dynamics and ageing-related diseases in humans. Philosophical Transactions of the Royal Society of London. Series B, Biological Sciences, 373(1741), 20160436. https://doi.org/10.1098/rstb.2016.0436
Baerlocher, G. M., Rice, K., Vulto, I., & Lansdorp, P. M. (2007). Longitudinal data on telomere length in leukocytes from newborn baboons support a marked drop in stem cell turnover around 1 year of age. Aging Cell, 6(1), 121–123. https://doi.org/10.1111/j.1474-9726.2006.00254.x
Baird, D. M., Rowson, J., Wynford-Thomas, D., & Kipling, D. (2003). Extensive allelic variation and ultrashort telomeres in senescent human cells. Nature Genetics, 33(2), 203–207. https://doi.org/10.1038/ng1084
Bakaysa, S. L., Mucci, L. a., Slagboom, P. E., Boomsma, D. I., Mcclearn, G. E., Johansson, B., & Pedersen, N. L. (2007). Telomere length predicts survival independent of genetic influences. Aging Cell, 6(6), 769–774. https://doi.org/10.1111/j.1474-9726.2007.00340.x
Barrett, E. L. B., Burke, T. a., Hammers, M., Komdeur, J., & Richardson, D. S. (2013). Telomere length and dynamics predict mortality in a wild longitudinal study. Molecular Ecology, 22(1), 249–259. https://doi.org/10.1111/mec.12110
Bates, D., Mächler, M., Bolker, B., & Walker, S. (2015). Fitting Linear Mixed-Effects Models Using lme4. Journal of Statistical Software, 67(1). https://doi.org/10.18637/jss.v067.i01
Bateson, M, & Nettle, D. (2016). The telomere lengthening conundrum - it could be biology. Aging Cell, 2016, 1–8. https://doi.org/10.1111/acel.12555
Bateson, Melissa. (2016). Cumulative stress in research animals: Telomere attrition as a biomarker in a welfare context? BioEssays, 38(2), 201–212. https://doi.org/10.1002/bies.201500127
Benetos, A., Kark, J. D., Susser, E., Kimura, M., Sinnreich, R., Chen, W., … Aviv, A. (2013). Tracking and fixed ranking of leukocyte telomere length across the adult life course. Aging Cell, 12(4), 615–621. https://doi.org/10.1111/acel.12086
Bichet, C., Bouwhuis, S., Bauch, C., Verhulst, S., Becker, P. H., & Vedder, O. (2019). Telomere length is repeatable, shortens with age and reproductive success, and predicts remaining lifespan in a long‐lived seabird. Molecular Ecology, 2020(29), 429–441. https://doi.org/https://doi.org/10.1111/mec.1533
Bize, P., Criscuolo, F., Metcalfe, N. B., Nasir, L., & Monaghan, P. (2009). Telomere dynamics rather than age predict life expectancy in the wild. Proceedings of the Royal Society of London B, 276(1662), 1679–1683. https://doi.org/10.1098/rspb.2008.1817
Blackburn, E., & Gall, J. (1978). A Tandemly Repeated Sequence at the Termini of Extrachromosomal Ribbosomal RNA Genes in Tetrahymena. Journal of Molecular Biology, (120), 33–53.
Boonekamp, J., Mulder, G., Salomons, H. M., Dijkstra, C., & Verhulst, S. (2014). Nestling telomere shortening, but not telomere length, reflects developmental stress and predicts survival in wild birds. Proceedings of the Royal Society B: Biological Sciences, 281(April), 20133287. https://doi.org/http://dx.doi.org/10.1098/rspb.2013.3287
Boonekamp, J., Simons, M. J. P., Hemerik, L., & Verhulst, S. (2013). Telomere length behaves as biomarker of somatic redundancy rather than biological age. Aging Cell, 12(2), 330–332. https://doi.org/10.1111/acel.12050
Boonekamp, J., Mulder, G., Salomons, H., Dijkstra, C., & Verhulst, S. (2014). Nestling telomere shortening, but not telomere length, reflects developmental stress and predicts survival in wild birds. Proceedings of the Royal Society B: Biological Sciences, 281(April), 20133287. https://doi.org/http://dx.doi.org/10.1098/rspb.2013.3287
Broer, L., Codd, V., Nyholt, D. R., Deelen, J., Mangino, M., Willemsen, G., … Boomsma, D. I. (2013). Meta-analysis of telomere length in 19 713 subjects reveals high heritability, stronger maternal inheritance and a paternal age effect. European Journal of Human Genetics, 21(10), 1163–1168. https://doi.org/10.1038/ejhg.2012.303
Brown, D. E., Dechow, C. D., Liu, W. S., Harvatine, K. J., & Ott, T. L. (2012). Hot topic: Association of telomere length with age, herd, and culling in lactating Holsteins. Journal of Dairy Science, 95(11), 6384–6387. https://doi.org/10.3168/jds.2012-5593
Cawthon, R., Smith, K., O’Brien, E., Sivatchenko, A., & Kerber, R. (2003). Association between telomere length in blood and mortality in people aged 60 years or older. The Lancet, 361(9355), 393–395. https://doi.org/10.1016/S0140-6736(03)12384-7
Chen, W., Kimura, M., Kim, S., Cao, X., Srinivasan, S. R., Berenson, G. S., … Aviv, A. (2011). Longitudinal versus cross-sectional evaluations of leukocyte telomere length dynamics: Age-dependent telomere shortening is the rule. Journals of Gerontology - Series A Biological Sciences and Medical Sciences, 66 A(3), 312–319. https://doi.org/10.1093/gerona/glq223
Cherkas, L F, Aviv, a, Valdes, a M., Hunkin, J. L., Gardner, J. P., Surdulescu, G. L., … Spector, T. D. (2006). The effects of social status on biological aging as measured by white-blood-cell telomere length. Aging Cell, 5(5), 361–365. https://doi.org/10.1111/j.1474-9726.2006.00222.x
Cherkas, Lynn F, Hunkin, J. L., Kato, B. S., Richards, J. B., Gardner, J. P., Surdulescu, G. L., … Aviv, A. (2008). The association between physical activity in leisure time and leukocyte telomere length. Archives of Internal Medicine, 168(2), 154–158. https://doi.org/10.1001/archinternmed.2007.39
De Lange, T. (2005). Shelterin: The protein complex that shapes and safeguards human telomeres. Genes and Development, 19(18), 2100–2110. https://doi.org/10.1101/gad.1346005
Dugdale, H. L., & Richardson, D. S. (2018). Heritability of telomere variation: it’s all about the environment! Philosophical Transactions of the Royal Society B, 373(1741), 20160450. https://doi.org/10.1098/rstb.2016.0450
Entringer, S., Epel, E. S., Kumsta, R., Lin, J., Hellhammer, D. H., Blackburn, E. H., … Wadhwa, P. D. (2011). Stress exposure in intrauterine life is associated with shorter telomere length in young adulthood. Proceedings of the National Academy of Sciences, 108(33), E513–E518. https://doi.org/10.1073/pnas.1107759108
Entringer, Sonja, Epel, E. S., Lin, J., Buss, C., Shahbaba, B., Blackburn, E. H., … Wadhwa, P. D. (2013). Maternal psychosocial stress during pregnancy is associated with newborn leukocyte telomere length. American Journal of Obstetrics and Gynecology, 208(2), 134.e1-134.e7. https://doi.org/10.1016/j.ajog.2012.11.033
Epel, E., Daubenmier, J., Moskowitz, J. T., Folkman, S., & Blackburn, E. (2009). Can meditation slow rate of cellular aging? Cognitive stress, mindfulness, and telomeres. Annals of the New York Academy of Sciences, 1172, 34–53. https://doi.org/10.1111/j.1749-6632.2009.04414.x
Epel, E. S., Blackburn, E. H., Lin, J., Dhabhar, F. S., Adler, N. E., Morrow, J. D., & Cawthon, R. M. (2004). Accelerated telomere shortening in response to life stress. Proceedings of the National Academy of Sciences, 101(49), 17312–17315. https://doi.org/10.1073/pnas.0407162101
Fairlie, J., Holland, R., Pilkington, J. G., Pemberton, J. M., Harrington, L., & Nussey, D. H. (2015). Lifelong leukocyte telomere dynamics and survival in a free-living mammal. Aging Cell, 140–148. https://doi.org/10.1111/acel.12417
Farzaneh-Far, R., Lin, J., Epel, E., Lapham, K., Blackburn, E., & Whooley, M. A. (2010). Telomere length trajectory and its determinants in persons with coronary artery disease: Longitudinal findings from the heart and soul study. PLoS ONE, 5(1). https://doi.org/10.1371/journal.pone.0008612
Gardner, J. P., Li, S., Srinivasan, S. R., Chen, W., Kimura, M., Lu, X., … Aviv, A. (2005). Rise in insulin resistance is associated with escalated telomere attrition. Circulation, 111(17), 2171–2177. https://doi.org/10.1161/01.CIR.0000163550.70487.0B
Goglin, S. E., Farzaneh-Far, R., Epel, E. S., Lin, J., Blackburn, E. H., & Whooley, M. A. (2016). Change in leukocyte telomere length predicts mortality in patients with stable coronary heart disease from the heart and soul study. PLoS ONE, 11(10), 1–12. https://doi.org/10.1371/journal.pone.0160748
Gomes, N. M. V, Ryder, O. a., Houck, M. L., Charter, S. J., Walker, W., Forsyth, N. R., … Wright, W. E. (2011). Comparative biology of mammalian telomeres: Hypotheses on ancestral states and the roles of telomeres in longevity determination. Aging Cell, 10(5), 761–768. https://doi.org/10.1111/j.1474-9726.2011.00718.x
Harley, C., Futcher, A., & Greider, C. (1990). Telomeres shorten during ageing of human fibroblasts. Nature, 345(6274), 458–460. https://doi.org/10.1038/345458a0
Harrington, L., & Pucci, F. (2018). In medio stat virtus: unanticipated consequences of telomere dysequilibrium. Philosophical Transactions of the Royal Society B: Biological Sciences, 373(1741), 20160444. https://doi.org/10.1098/rstb.2016.0444
Haussmann, M. F., Longenecker, a. S., Marchetto, N. M., Juliano, S. a., & Bowden, R. M. (2012). Embryonic exposure to corticosterone modifies the juvenile stress response, oxidative stress and telomere length. Proceedings of the Royal Society B: Biological Sciences, 279(1732), 1447–1456. https://doi.org/10.1098/rspb.2011.1913
Haussmann, Mark F, Winkler, D. W., & Vleck, C. M. (2005). Longer telomeres associated with higher survival in birds. Biology Letters, 1(2), 212–214. https://doi.org/10.1098/rsbl.2005.0301
Heidinger, B. J., Blount, J. D., Boner, W., Griffiths, K., Metcalfe, N. B., & Monaghan, P. (2012). Telomere length in early life predicts lifespan. Proceedings of the National Academy of Sciences, 109(5), 1743–1748. https://doi.org/10.1073/pnas.1113306109
Huzen, J., Wong, L. S. M., van Veldhuisen, D. J., Samani, N. J., Zwinderman, A. H., Codd, V., … Van der Harst, P. (2014). Telomere length loss due to smoking and metabolic traits. Journal of Internal Medicine, 275(2), 155–163. https://doi.org/10.1111/joim.12149
Kark, J., Goldberger, N., Kimura, M., Sinnreich, R., & Aviv, A. (2012). Energy intake and leukocyte telomere length in young adults. The American Journal of Clinical Nutrition, (95), 479–487. https://doi.org/10.3945/ajcn.111.024521.1
Kim, S., Parks, C. G., DeRoo, L. A., Chen, H., Taylor, J. A., Cawthon, R. M., & Sandler, D. P. (2009). Obesity and weight gain in adulthood and telomere length. Cancer Epidemiology Biomarkers and Prevention, 18(3), 816–820. https://doi.org/10.1158/1055-9965.EPI-08-0935
Lai, T.-P., Zhang, N., Noh, J., Mender, I., Tedone, E., Huang, E., … Shay, J. W. (2017). A method for measuring the distribution of the shortest telomeres in cells and tissues. Nature Communications, 8(1), 1356. https://doi.org/10.1038/s41467-017-01291-z
Muñoz-lorente, M. A., Cano-martin, A. C., & Blasco, M. A. (2019). Mice with hyper-long telomeres show less metabolic aging and longer lifespans. Nature Communications, 10(4723), 1–14. https://doi.org/10.1038/s41467-019-12664-x
Nettle, D., Andrews, C., Reichert, S., Bedford, T., Kolenda, C., Parker, C., … Bateson, M. (2017). Early-life adversity accelerates cellular ageing and affects adult inflammation: Experimental evidence from the European starling. Scientific Reports, 7: 40794(December 2016), 1–10. https://doi.org/10.1038/srep40794
Nettle, D., Monaghan, P., Boner, W., Gillespie, R., & Bateson, M. (2013). Bottom of the heap: Having heavier competitors accelerates early-life telomere loss in the European starling, Sturnus vulgaris. PLoS ONE, 8(12). https://doi.org/10.1371/journal.pone.0083617
Nettle, D., Seeker, L., Nussey, D., Froy, H., & Bateson, M. (2019). Consequences of measurement error in qPCR telomere data: A simulation study. PLOS ONE, 14(5), e0216118. https://doi.org/10.1371/journal.pone.0216118
Nordfjäll, K., Svenson, U., Norrback, K. F., Adolfsson, R., Lenner, P., & Roos, G. (2009). The individual blood cell telomere attrition rate is telomere length dependent. PLoS Genetics, 5(2), 2–7. https://doi.org/10.1371/journal.pgen.1000375
Olovnikov, A. M. (1973). A theory of marginotomy. The incomplete copying of template margin in enzymic synthesis of polynucleotides and biological significance of the phenomenon. Journal of Theoretical Biology, 41(1), 181–190. https://doi.org/10.1016/0022-5193(73)90198-7
R Core Team. (2014). R: A Language and Environment for Statistical Computing. Vienna, Austria: R Foundation for Statistical Computing. Retrieved from http://www.r-project.org/
Reichert, S., Stier, A., Zahn, S., Arrivé, M., Bize, P., Massemin, S., & Criscuolo, F. (2014). Increased brood size leads to persistent eroded telomeres. Frontiers in Ecology and Evolution, 2(April), 1–11. https://doi.org/10.3389/fevo.2014.00009
Seeker, L. A., Holland, R., Underwood, S., Fairlie, J., Psifidi, A., Ilska, J., … Nussey, D. H. (2016). Method specific calibration corrects for DNA extraction method effects on relative telomere length measurements by quantitative PCR. PLoS ONE, 11(10), 1–15. https://doi.org/10.1371/journal.pone.0164046
Seeker, L. A., Ilska, J., Psifidi, A., Wilbourn, R. V., Underwood, S. L., Fairlie, J., … Banos, G. (2018). Bovine telomere dynamics and the association between telomere length and productive lifespan. Scientific Reports.
Seeker, L. A., Ilska, J., Psifidi, A., Wilbourn, R. V, Underwood, L., Fairlie, J., … Banos, G. (2018). Longitudinal changes in telomere length and associated genetic parameters in dairy cattle analysed using random regression models. Plos One, 1–15. https://doi.org/10.1371/journal.pone.0192864
Seeker, L., Holland, B., Psifidi, A., Banos, B., & Nussey, D. (2015). A robust assay for the measurement of telomere length in dairy cattle. In BSAS annual conference (p. 066).
Shalev, I., Moffitt, T. E., Sugden, K., Williams, B., Houts, R. M., Danese, A., … Caspi, A. (2012). Exposure to violence during childhood is associated with telomere erosion from 5 to 10 years of age: a longitudinal study. Mol.Psychiatry, 18(5), 576–581. https://doi.org/10.1038/mp.2012.32
Spurgin, L. G., Bebbington, K., Fairfield, E. A., Hammers, M., Komdeur, J., Burke, T., … Richardson, D. S. (2018). Spatio-temporal variation in lifelong telomere dynamics in a long-term ecological study. Journal of Animal Ecology, 87(1), 187–198. https://doi.org/10.1111/1365-2656.12741
Steenstrup, T., Hjelmborg, J. V. B., Mortensen, L. H., Kimura, M., Christensen, K., & Aviv, A. (2013). Leukocyte telomere dynamics in the elderly. European Journal of Epidemiology, 28(2), 181–187. https://doi.org/10.1007/s10654-013-9780-4
Svenson, U., Nordfjäll, K., Baird, D., Roger, L., Osterman, P., Hellenius, M.-L., & Roos, G. (2011). Blood Cell Telomere Length Is a Dynamic Feature. PLOS ONE, 6(6), e21485. Retrieved from https://doi.org/10.1371/journal.pone.0021485
Therneau, T. M. (2015). A Package for Survival Analysis in S. Retrieved from https://cran.r-project.org/package=survival
Veerkamp, R. ., Simm, G., & Oldham, J. . (1994). Effects of interaction between genotype and feeding system on milk production, feed intake, efficiency and body tissue mobilization in dairy cows. Livestock Production Science, 39(3), 229–241. https://doi.org/10.1016/0301-6226(94)90202-X
Voillemot, M., Hine, K., Zahn, S., Criscuolo, F., Gustafsson, L., Doligez, B., & Bize, P. (2012). Effects of brood size manipulation and common origin on phenotype and telomere length in nestling collared flycatchers. BMC Ecology, 12(1), 17. https://doi.org/10.1186/1472-6785-12-17
Watson, J. D. (1972). Origin of concatemeric T7 DNA. Nature: New Biology, 239(94), 197–201. Retrieved from http://www.ncbi.nlm.nih.gov/pubmed/4507727
Whittemore, K., Vera, E., Martínez-nevado, E., Sanpera, C., & Blasco, M. A. (2019). Telomere shortening rate predicts species life span. Proceedings of the National Academy of Sciences, 116(30), 15122–15127. https://doi.org/10.1073/pnas.1902452116
Wickham, H. (2009). ggplot2: Elegant Graphics for Data Analysis. Springer-Verlag New York. Retrieved from http://ggplot2.org
Wilbourn, R. V., Moatt, J. P., Froy, H., Walling, C. A., Nussey, D. H., & Boonekamp, J. J. (2018). The relationship between telomere length and mortality risk in non-model vertebrate systems: a meta-analysis. Phil. Trans. R. Soc. B, 373(1741), 20160447. https://doi.org/10.1098/RSTB.2016.0447
Yuan, X., Kronström, M., Hellenius, M. L., Cederholm, T., & Xu, D. (2018). Longitudinal changes in leukocyte telomere length and mortality in elderly Swedish men. Aging, 10(10), 3005–3016.
Data Accessibility
The complete dataset will be accessible after acceptance on GitHub.
PART 2: THE ANALYSIS
Load libraries
library(ggplot2)
library(nlme)
library(lme4)## Loading required package: Matrix##
## Attaching package: 'lme4'## The following object is masked from 'package:nlme':
##
## lmListlibrary(lmerTest)##
## Attaching package: 'lmerTest'## The following object is masked from 'package:lme4':
##
## lmer## The following object is masked from 'package:stats':
##
## steplibrary(grid)
library(gridExtra)
library(mgcv) # for GAM## This is mgcv 1.8-31. For overview type 'help("mgcv-package")'.library(data.table)
library(stringr)
library(magrittr)
library(ggthemes)
library(plyr)
library(dplyr)##
## Attaching package: 'dplyr'## The following objects are masked from 'package:plyr':
##
## arrange, count, desc, failwith, id, mutate, rename, summarise,
## summarize## The following objects are masked from 'package:data.table':
##
## between, first, last## The following object is masked from 'package:gridExtra':
##
## combine## The following object is masked from 'package:nlme':
##
## collapse## The following objects are masked from 'package:stats':
##
## filter, lag## The following objects are masked from 'package:base':
##
## intersect, setdiff, setequal, unionlibrary(survival)
library(survminer)## Loading required package: ggpubr## Registered S3 methods overwritten by 'car':
## method from
## influence.merMod lme4
## cooks.distance.influence.merMod lme4
## dfbeta.influence.merMod lme4
## dfbetas.influence.merMod lme4##
## Attaching package: 'ggpubr'## The following object is masked from 'package:plyr':
##
## mutatelibrary(ggsci)
#library(AICcmodavg)Build colour palettes
# Make colour palette
mypal <- pal_npg("nrc", alpha = 0.7)(10)
mypal2 <-pal_tron("legacy", alpha = 0.7)(7)
mypal3<-pal_lancet("lanonc", alpha = 0.7)(9)
mypal4<-pal_simpsons(palette = c("springfield"), alpha = 0.7)(16)
mycoloursP<- c(mypal, mypal2, mypal3, mypal4)Read in data
data<- read.csv("/Users/luise/Documents/Work/CattleTelomereLength/LongTelomereChangePaper/git/TelomereChangeInDairyCattle/LongTelomereDataCattle.csv")
# create data which only lists each animal ID once (needed for some counts/ histograms below)
uData<- data[!duplicated(data$recoded_id),]
# create data that excludes all first measurements without t-1 time point. This is used for investigations of change
redData<-subset(data, data$residual_rltl_tminus1 != "NULL")
init_change_dat<-read.csv("/Users/luise/Documents/Work/CattleTelomereLength/LongTelomereChangePaper/git/TelomereChangeInDairyCattle/first_year_long_tel_cattle.csv")
coxData<-read.csv("/Users/luise/Documents/Work/CattleTelomereLength/LongTelomereChangePaper/git/TelomereChangeInDairyCattle/cox_data.csv")
coxDataR<-subset(coxData, coxData$amp_305_lactation != "NULL")Fig. 1
cor.test(redData$residual_rltl_t,as.numeric(paste(redData$residual_rltl_tminus1)))##
## Pearson's product-moment correlation
##
## data: redData$residual_rltl_t and as.numeric(paste(redData$residual_rltl_tminus1))
## t = 13.058, df = 1018, p-value < 2.2e-16
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
## 0.3249385 0.4301473
## sample estimates:
## cor
## 0.3787659SubsPlot<-ggplot(redData, aes(x=as.numeric(paste(residual_rltl_tminus1)), y=redData$residual_rltl_t))+geom_point()+xlab("RLTL at time t-1")+ylab("RLTL at time t")+theme_classic(20)+geom_abline(intercept=0, slope=1, colour="red")+ stat_smooth(method=lm, formula=y~x)+
annotate("text", x = -0.4, y = 0.5, label = "(B)", fontface =2, size = 10)
agePlot<- ggplot(redData, aes(x=as.factor(redData$age_y), y=as.numeric(paste(redData$diff_rltl_res)), colour = as.factor(redData$age_y)))+geom_jitter()+theme_classic(20)+ xlab("Age in years")+ ylab("RLTL change")+ geom_boxplot()+
annotate("text", x = 1, y = 0.5, label = "(A)", fontface =2, size = 10) + scale_color_manual(values=mycoloursP[16:40])+ theme(legend.position="none")
grid.arrange(agePlot, SubsPlot, ncol=2)## Warning: Use of `redData$age_y` is discouraged. Use `age_y` instead.## Warning: Use of `redData$diff_rltl_res` is discouraged. Use `diff_rltl_res`
## instead.## Warning: Use of `redData$age_y` is discouraged. Use `age_y` instead.
## Warning: Use of `redData$age_y` is discouraged. Use `age_y` instead.## Warning: Use of `redData$diff_rltl_res` is discouraged. Use `diff_rltl_res`
## instead.## Warning: Use of `redData$age_y` is discouraged. Use `age_y` instead.## Warning: Use of `redData$residual_rltl_t` is discouraged. Use `residual_rltl_t`
## instead.
## Warning: Use of `redData$residual_rltl_t` is discouraged. Use `residual_rltl_t`
## instead.
Fig. 2
surv_plotlist<-list()
surv_plotlist[[1]] <- init_ch_plot
surv_plotlist[[2]] <- mean_abs_rltl_change
surv_plotlist[[3]] <- meanRLTL
surv_plotlist[[4]] <- mean_rltl_change
arrange_ggsurvplots(surv_plotlist, print = TRUE, ncol = 2, nrow = 2)
calculate sample number per animal
count<- nrow(data)
unique_count<- length(unique(data$recoded_id))
sample_year_summ<-summary(as.factor(data$sample_year))
print(paste("There are ", count, " relative leukocyte telomere length (RLTL) measurements of ", unique_count, " animals included in the study.", sep =""))## [1] "There are 1325 relative leukocyte telomere length (RLTL) measurements of 305 animals included in the study."Distibution of sample years
sample_year_summ## 2008 2009 2010 2011 2012 2013 2014
## 22 91 179 247 266 272 248Plot sample number per animal
hist1<-ggplot(uData, aes(x = uData$SamplesPerAnimal)) +
geom_histogram(colour= "black", fill = mycoloursP[7], binwidth = 1) + xlab("Number of samples per animal") +
geom_vline(xintercept= mean(as.numeric(paste(uData$SamplesPerAnimal))), colour="red")+
annotate("text", x = 1, y = 100, label = "(A)", fontface =2, size = 8) + theme_classic(20)+
scale_x_continuous(breaks=c(1:8))
#hist1Plot distribution of telomere residuals
hist2<-ggplot(data, aes(x = data$residual_rltl_t)) +
geom_histogram(colour= "black", fill = mycoloursP[12]) + xlab("Adjusted RLTL") +
geom_vline(xintercept= mean(as.numeric(paste(data$residual_rltl_t))), colour="red") +
annotate("text", x = -0.4, y = 120, label = "(B)", fontface =2, size = 8)+
theme_classic(20)+
stat_function(
fun = function(x, mean, sd, n){
n/25* dnorm(x = x, mean = mean, sd = sd)
},
args = with(data, c(mean = mean(data$residual_rltl_t), sd = sd(data$residual_rltl_t), n
= length(data$residual_rltl_t)))
)
#hist2# statistical tests for normal distribution of RTL measurements
shapiro.test(data$residual_rltl_t)##
## Shapiro-Wilk normality test
##
## data: data$residual_rltl_t
## W = 0.989, p-value = 1.967e-08ks.test(data$residual_rltl_t, pnorm)##
## One-sample Kolmogorov-Smirnov test
##
## data: data$residual_rltl_t
## D = 0.36744, p-value < 2.2e-16
## alternative hypothesis: two-sidedqqnorm(data$residual_rltl_t)
Plot difference between RLTL res at time t and time t-1
diff<- as.numeric(paste(redData$diff_rltl_res))
hist3<-ggplot(redData, aes(x = as.numeric(paste(redData$diff_rltl_res)))) +
geom_histogram(colour= "black", fill = mycoloursP[13]) + xlab("Adjusted RLTL change") +
geom_vline(xintercept= mean(as.numeric(paste(redData$diff_rltl_res))), colour="red") +
annotate("text", x = -0.6, y = 95, label = "(C)", fontface =2, size= 8)+
theme_classic(20)+
stat_function(
fun = function(x, mean, sd, n){
n/23* dnorm(x = x, mean = mean, sd = sd)
},
args = with(redData, c(mean = mean(as.numeric(paste(redData$diff_rltl_res))), sd = sd(as.numeric(paste(redData$diff_rltl_res))), n
= length(as.numeric(paste(redData$diff_rltl_res))))))
#hist3grid.arrange(hist1, hist2, hist3, ncol = 3)## Warning: Use of `uData$SamplesPerAnimal` is discouraged. Use `SamplesPerAnimal`
## instead.## Warning: Use of `data$residual_rltl_t` is discouraged. Use `residual_rltl_t`
## instead.
## Warning: Use of `data$residual_rltl_t` is discouraged. Use `residual_rltl_t`
## instead.## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.## Warning: Use of `redData$diff_rltl_res` is discouraged. Use `diff_rltl_res`
## instead.## Warning: Use of `redData$diff_rltl_res` is discouraged. Use `diff_rltl_res`
## instead.## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
# statistical tests for normal distribution of difference in RTL residuals
shapiro.test(diff)##
## Shapiro-Wilk normality test
##
## data: diff
## W = 0.998, p-value = 0.2655ks.test(diff, pnorm)##
## One-sample Kolmogorov-Smirnov test
##
## data: diff
## D = 0.35414, p-value < 2.2e-16
## alternative hypothesis: two-sidedqqnorm(diff)
# test if mean RLTL change is significantly different from zero:
t.test(as.numeric(paste(redData$diff_rltl_res)), mu=0, alternative="less", conf.level=0.99)##
## One Sample t-test
##
## data: as.numeric(paste(redData$diff_rltl_res))
## t = -5.8438, df = 1019, p-value = 3.43e-09
## alternative hypothesis: true mean is less than 0
## 99 percent confidence interval:
## -Inf -0.01949685
## sample estimates:
## mean of x
## -0.03242526Shortening and elongating:
shortening<-subset(redData, as.numeric(paste(redData$diff_rltl_res)) < 0)
elongating<-subset(redData, as.numeric(paste(redData$diff_rltl_res))>0)
percent_shortening<-nrow(shortening)/nrow(redData)*100
percent_shortening## [1] 56.47059percent_elongating<-nrow(elongating)/nrow(redData)*100
percent_elongating## [1] 43.52941Initial full model of RLTL change
init_mod<-lmer(as.numeric(paste(redData$diff_rltl_res))~ redData$age_y + redData$feed_group + redData$genetic_group+ as.factor(redData$birth_year)+ as.factor(redData$sample_year) + as.numeric(paste(redData$sample_interval))+ redData$health_event_2week +(1|redData$recoded_id), na.action=na.exclude)## boundary (singular) fit: see ?isSingularsummary(init_mod)## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula:
## as.numeric(paste(redData$diff_rltl_res)) ~ redData$age_y + redData$feed_group +
## redData$genetic_group + as.factor(redData$birth_year) + as.factor(redData$sample_year) +
## as.numeric(paste(redData$sample_interval)) + redData$health_event_2week +
## (1 | redData$recoded_id)
##
## REML criterion at convergence: -599.3
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -3.4980 -0.6827 0.0248 0.6712 3.5498
##
## Random effects:
## Groups Name Variance Std.Dev.
## redData$recoded_id (Intercept) 0.00000 0.0000
## Residual 0.02953 0.1718
## Number of obs: 1020, groups: redData$recoded_id, 305
##
## Fixed effects:
## Estimate Std. Error df
## (Intercept) 2.376e-01 1.737e-01 1.002e+03
## redData$age_y 2.219e-02 1.385e-02 1.002e+03
## redData$feed_group2 -8.381e-03 1.223e-02 1.002e+03
## redData$feed_groupNULL -6.651e-03 1.627e-02 1.002e+03
## redData$genetic_groupS -1.025e-02 1.144e-02 1.002e+03
## as.factor(redData$birth_year)2009 -1.431e-02 2.062e-02 1.002e+03
## as.factor(redData$birth_year)2010 -1.231e-02 2.855e-02 1.002e+03
## as.factor(redData$birth_year)2011 -6.343e-03 3.863e-02 1.002e+03
## as.factor(redData$birth_year)2012 -2.192e-02 5.533e-02 1.002e+03
## as.factor(redData$birth_year)2013 -1.466e-01 1.114e-01 1.002e+03
## as.factor(redData$birth_year)2014 -2.415e-01 1.925e-01 1.002e+03
## as.factor(redData$sample_year)2010 -3.639e-01 1.748e-01 1.002e+03
## as.factor(redData$sample_year)2011 -3.223e-01 1.763e-01 1.002e+03
## as.factor(redData$sample_year)2012 -3.008e-01 1.784e-01 1.002e+03
## as.factor(redData$sample_year)2013 -2.767e-01 1.824e-01 1.002e+03
## as.factor(redData$sample_year)2014 -3.003e-01 1.886e-01 1.002e+03
## as.numeric(paste(redData$sample_interval)) -2.685e-06 5.434e-05 1.002e+03
## redData$health_event_2week 1.764e-02 2.447e-02 1.002e+03
## t value Pr(>|t|)
## (Intercept) 1.368 0.1715
## redData$age_y 1.602 0.1095
## redData$feed_group2 -0.685 0.4933
## redData$feed_groupNULL -0.409 0.6828
## redData$genetic_groupS -0.896 0.3706
## as.factor(redData$birth_year)2009 -0.694 0.4880
## as.factor(redData$birth_year)2010 -0.431 0.6663
## as.factor(redData$birth_year)2011 -0.164 0.8696
## as.factor(redData$birth_year)2012 -0.396 0.6920
## as.factor(redData$birth_year)2013 -1.316 0.1886
## as.factor(redData$birth_year)2014 -1.255 0.2099
## as.factor(redData$sample_year)2010 -2.082 0.0376 *
## as.factor(redData$sample_year)2011 -1.829 0.0677 .
## as.factor(redData$sample_year)2012 -1.686 0.0921 .
## as.factor(redData$sample_year)2013 -1.517 0.1295
## as.factor(redData$sample_year)2014 -1.592 0.1117
## as.numeric(paste(redData$sample_interval)) -0.049 0.9606
## redData$health_event_2week 0.721 0.4711
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1##
## Correlation matrix not shown by default, as p = 18 > 12.
## Use print(x, correlation=TRUE) or
## vcov(x) if you need it## convergence code: 0
## boundary (singular) fit: see ?isSingularanova(init_mod)## Type III Analysis of Variance Table with Satterthwaite's method
## Sum Sq Mean Sq NumDF DenDF F value
## redData$age_y 0.07576 0.075762 1 1002 2.5657
## redData$feed_group 0.01465 0.007327 2 1002 0.2481
## redData$genetic_group 0.02369 0.023694 1 1002 0.8024
## as.factor(redData$birth_year) 0.12877 0.021462 6 1002 0.7268
## as.factor(redData$sample_year) 0.42749 0.085498 5 1002 2.8954
## as.numeric(paste(redData$sample_interval)) 0.00007 0.000072 1 1002 0.0024
## redData$health_event_2week 0.01535 0.015351 1 1002 0.5199
## Pr(>F)
## redData$age_y 0.10952
## redData$feed_group 0.78031
## redData$genetic_group 0.37059
## as.factor(redData$birth_year) 0.62806
## as.factor(redData$sample_year) 0.01329 *
## as.numeric(paste(redData$sample_interval)) 0.96060
## redData$health_event_2week 0.47107
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1AIC(init_mod)## [1] -559.2719# There is a warning that this model is singular which is due to the random effect variance structure. It means that the variance due to the animal ID is not significantly different from 0. Although this may be a problem for some models, I think in this case it is not. It simply means that there is no variation between animals in the RLTL change and that is to be expected here, because the variance within animals is much larger than the variance between animals.
# Sample year is statistically significant, but it might only be signifgicant due to the design of our study. Becasue we wanted to have repeat measurements for all animals with the first measurement taken with the first 15 days of life, sampling in the first year only included calves which have been shown before to have longer telomeres than adults and to shorten them in their first year of life. So is the year effect simply driven by that? If so, it should vanish after excluding samples taken within the two yearsFull model of RLTL change excluding first two years
summary(as.factor(data$sample_year))## 2008 2009 2010 2011 2012 2013 2014
## 22 91 179 247 266 272 248test_dat<-subset(data, data$sample_year != 2008)
test_dat<-subset(test_dat, test_dat$sample_year != 2009)
long_dat<-subset(test_dat, test_dat$diff_rltl_res != "NULL")
late_mod<-lmer(as.numeric(paste(long_dat$diff_rltl_res))~ long_dat$age_y + long_dat$feed_group + long_dat$genetic_group + as.factor(long_dat$birth_year)+ as.factor(long_dat$sample_year) + as.numeric(paste(long_dat$sample_interval))+ long_dat$health_event_2week +(1|long_dat$recoded_id), na.action=na.exclude)## boundary (singular) fit: see ?isSingularsummary(late_mod)## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: as.numeric(paste(long_dat$diff_rltl_res)) ~ long_dat$age_y +
## long_dat$feed_group + long_dat$genetic_group + as.factor(long_dat$birth_year) +
## as.factor(long_dat$sample_year) + as.numeric(paste(long_dat$sample_interval)) +
## long_dat$health_event_2week + (1 | long_dat$recoded_id)
##
## REML criterion at convergence: -599.3
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -3.4980 -0.6836 0.0252 0.6717 3.5498
##
## Random effects:
## Groups Name Variance Std.Dev.
## long_dat$recoded_id (Intercept) 0.00000 0.0000
## Residual 0.02953 0.1718
## Number of obs: 1019, groups: long_dat$recoded_id, 305
##
## Fixed effects:
## Estimate Std. Error df
## (Intercept) -1.263e-01 3.879e-02 1.002e+03
## long_dat$age_y 2.219e-02 1.385e-02 1.002e+03
## long_dat$feed_group2 -8.381e-03 1.223e-02 1.002e+03
## long_dat$feed_groupNULL -6.651e-03 1.627e-02 1.002e+03
## long_dat$genetic_groupS -1.025e-02 1.144e-02 1.002e+03
## as.factor(long_dat$birth_year)2009 -1.431e-02 2.062e-02 1.002e+03
## as.factor(long_dat$birth_year)2010 -1.231e-02 2.855e-02 1.002e+03
## as.factor(long_dat$birth_year)2011 -6.343e-03 3.863e-02 1.002e+03
## as.factor(long_dat$birth_year)2012 -2.192e-02 5.533e-02 1.002e+03
## as.factor(long_dat$birth_year)2013 -1.466e-01 1.114e-01 1.002e+03
## as.factor(long_dat$birth_year)2014 -2.415e-01 1.925e-01 1.002e+03
## as.factor(long_dat$sample_year)2011 4.158e-02 2.644e-02 1.002e+03
## as.factor(long_dat$sample_year)2012 6.309e-02 3.329e-02 1.002e+03
## as.factor(long_dat$sample_year)2013 8.721e-02 4.463e-02 1.002e+03
## as.factor(long_dat$sample_year)2014 6.360e-02 6.047e-02 1.002e+03
## as.numeric(paste(long_dat$sample_interval)) -2.685e-06 5.434e-05 1.002e+03
## long_dat$health_event_2week 1.764e-02 2.447e-02 1.002e+03
## t value Pr(>|t|)
## (Intercept) -3.256 0.00117 **
## long_dat$age_y 1.602 0.10952
## long_dat$feed_group2 -0.685 0.49328
## long_dat$feed_groupNULL -0.409 0.68275
## long_dat$genetic_groupS -0.896 0.37059
## as.factor(long_dat$birth_year)2009 -0.694 0.48796
## as.factor(long_dat$birth_year)2010 -0.431 0.66631
## as.factor(long_dat$birth_year)2011 -0.164 0.86961
## as.factor(long_dat$birth_year)2012 -0.396 0.69205
## as.factor(long_dat$birth_year)2013 -1.316 0.18862
## as.factor(long_dat$birth_year)2014 -1.255 0.20990
## as.factor(long_dat$sample_year)2011 1.573 0.11605
## as.factor(long_dat$sample_year)2012 1.895 0.05835 .
## as.factor(long_dat$sample_year)2013 1.954 0.05097 .
## as.factor(long_dat$sample_year)2014 1.052 0.29309
## as.numeric(paste(long_dat$sample_interval)) -0.049 0.96060
## long_dat$health_event_2week 0.721 0.47107
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1##
## Correlation matrix not shown by default, as p = 17 > 12.
## Use print(x, correlation=TRUE) or
## vcov(x) if you need it## convergence code: 0
## boundary (singular) fit: see ?isSingularanova(late_mod)## Type III Analysis of Variance Table with Satterthwaite's method
## Sum Sq Mean Sq NumDF DenDF
## long_dat$age_y 0.075762 0.075762 1 1002
## long_dat$feed_group 0.014654 0.007327 2 1002
## long_dat$genetic_group 0.023694 0.023694 1 1002
## as.factor(long_dat$birth_year) 0.128774 0.021462 6 1002
## as.factor(long_dat$sample_year) 0.312340 0.078085 4 1002
## as.numeric(paste(long_dat$sample_interval)) 0.000072 0.000072 1 1002
## long_dat$health_event_2week 0.015351 0.015351 1 1002
## F value Pr(>F)
## long_dat$age_y 2.5657 0.10952
## long_dat$feed_group 0.2481 0.78031
## long_dat$genetic_group 0.8024 0.37059
## as.factor(long_dat$birth_year) 0.7268 0.62806
## as.factor(long_dat$sample_year) 2.6443 0.03235 *
## as.numeric(paste(long_dat$sample_interval)) 0.0024 0.96060
## long_dat$health_event_2week 0.5199 0.47107
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1AIC(late_mod)## [1] -561.2719# sample year remains statistically significant. Reduced model of RLTL change
red_mod<-lmer(as.numeric(paste(redData$diff_rltl_res))~ redData$age_y + redData$feed_group + redData$genetic_group+ as.factor(redData$sample_year) + as.numeric(paste(redData$sample_interval)) + (1|redData$recoded_id), na.action=na.exclude)## boundary (singular) fit: see ?isSingularsummary(red_mod)## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula:
## as.numeric(paste(redData$diff_rltl_res)) ~ redData$age_y + redData$feed_group +
## redData$genetic_group + as.factor(redData$sample_year) +
## as.numeric(paste(redData$sample_interval)) + (1 | redData$recoded_id)
##
## REML criterion at convergence: -627.8
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -3.5005 -0.6953 0.0197 0.6760 3.5199
##
## Random effects:
## Groups Name Variance Std.Dev.
## redData$recoded_id (Intercept) 0.00000 0.0000
## Residual 0.02947 0.1717
## Number of obs: 1020, groups: redData$recoded_id, 305
##
## Fixed effects:
## Estimate Std. Error df
## (Intercept) 2.278e-01 1.724e-01 1.009e+03
## redData$age_y 2.571e-02 6.026e-03 1.009e+03
## redData$feed_group2 -1.045e-02 1.206e-02 1.009e+03
## redData$feed_groupNULL -1.499e-02 1.546e-02 1.009e+03
## redData$genetic_groupS -1.262e-02 1.134e-02 1.009e+03
## as.factor(redData$sample_year)2010 -3.693e-01 1.738e-01 1.009e+03
## as.factor(redData$sample_year)2011 -3.298e-01 1.736e-01 1.009e+03
## as.factor(redData$sample_year)2012 -3.087e-01 1.735e-01 1.009e+03
## as.factor(redData$sample_year)2013 -2.882e-01 1.740e-01 1.009e+03
## as.factor(redData$sample_year)2014 -3.190e-01 1.750e-01 1.009e+03
## as.numeric(paste(redData$sample_interval)) 1.029e-05 5.207e-05 1.009e+03
## t value Pr(>|t|)
## (Intercept) 1.321 0.1868
## redData$age_y 4.266 2.18e-05 ***
## redData$feed_group2 -0.867 0.3864
## redData$feed_groupNULL -0.969 0.3327
## redData$genetic_groupS -1.113 0.2660
## as.factor(redData$sample_year)2010 -2.125 0.0338 *
## as.factor(redData$sample_year)2011 -1.900 0.0578 .
## as.factor(redData$sample_year)2012 -1.779 0.0755 .
## as.factor(redData$sample_year)2013 -1.657 0.0979 .
## as.factor(redData$sample_year)2014 -1.823 0.0686 .
## as.numeric(paste(redData$sample_interval)) 0.198 0.8434
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) rdDt$_ rdD$_2 rD$_NU rdD$_S a.(D$_)2010 a.(D$_)2011
## redData$g_y 0.001
## rdDt$fd_gr2 -0.068 -0.099
## rdDt$f_NULL -0.042 0.114 0.370
## rdDt$gntc_S -0.064 0.090 -0.029 0.243
## a.(D$_)2010 -0.987 -0.034 0.032 0.010 0.021
## a.(D$_)2011 -0.988 -0.052 0.034 0.001 0.018 0.991
## a.(D$_)2012 -0.989 -0.067 0.037 0.006 0.019 0.992 0.995
## a.(D$_)2013 -0.986 -0.094 0.038 -0.001 0.015 0.991 0.995
## a.(D$_)2014 -0.980 -0.130 0.039 -0.005 0.009 0.988 0.992
## as.n((D$_)) -0.008 -0.033 0.047 -0.034 0.071 -0.099 -0.105
## a.(D$_)2012 a.(D$_)2013 a.(D$_)2014
## redData$g_y
## rdDt$fd_gr2
## rdDt$f_NULL
## rdDt$gntc_S
## a.(D$_)2010
## a.(D$_)2011
## a.(D$_)2012
## a.(D$_)2013 0.996
## a.(D$_)2014 0.994 0.995
## as.n((D$_)) -0.097 -0.105 -0.116
## convergence code: 0
## boundary (singular) fit: see ?isSingularanova(red_mod)## Type III Analysis of Variance Table with Satterthwaite's method
## Sum Sq Mean Sq NumDF DenDF F value
## redData$age_y 0.53630 0.53630 1 1009 18.2001
## redData$feed_group 0.03650 0.01825 2 1009 0.6193
## redData$genetic_group 0.03650 0.03650 1 1009 1.2388
## as.factor(redData$sample_year) 0.55730 0.11146 5 1009 3.7825
## as.numeric(paste(redData$sample_interval)) 0.00115 0.00115 1 1009 0.0390
## Pr(>F)
## redData$age_y 2.176e-05 ***
## redData$feed_group 0.538542
## redData$genetic_group 0.265975
## as.factor(redData$sample_year) 0.002125 **
## as.numeric(paste(redData$sample_interval)) 0.843439
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1AIC(red_mod) ## [1] -601.817Reduced model of RLTL change with age as two level factor
red_mod_fac<-lmer(as.numeric(paste(redData$diff_rltl_res))~ redData$AgeGroup + redData$feed_group + redData$genetic_group+ as.factor(redData$sample_year) + as.numeric(paste(redData$sample_interval))+(1|redData$recoded_id), na.action=na.exclude)## boundary (singular) fit: see ?isSingularsummary(red_mod_fac)## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: as.numeric(paste(redData$diff_rltl_res)) ~ redData$AgeGroup +
## redData$feed_group + redData$genetic_group + as.factor(redData$sample_year) +
## as.numeric(paste(redData$sample_interval)) + (1 | redData$recoded_id)
##
## REML criterion at convergence: -627
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -3.6485 -0.6861 0.0272 0.6780 3.3620
##
## Random effects:
## Groups Name Variance Std.Dev.
## redData$recoded_id (Intercept) 0.00000 0.0000
## Residual 0.02956 0.1719
## Number of obs: 1020, groups: redData$recoded_id, 305
##
## Fixed effects:
## Estimate Std. Error df
## (Intercept) 2.988e-01 1.737e-01 1.009e+03
## redData$AgeGroupyoung -7.262e-02 1.872e-02 1.009e+03
## redData$feed_group2 -5.832e-03 1.202e-02 1.009e+03
## redData$feed_groupNULL -1.786e-02 1.543e-02 1.009e+03
## redData$genetic_groupS -1.565e-02 1.131e-02 1.009e+03
## as.factor(redData$sample_year)2010 -3.503e-01 1.740e-01 1.009e+03
## as.factor(redData$sample_year)2011 -3.135e-01 1.737e-01 1.009e+03
## as.factor(redData$sample_year)2012 -2.886e-01 1.735e-01 1.009e+03
## as.factor(redData$sample_year)2013 -2.549e-01 1.737e-01 1.009e+03
## as.factor(redData$sample_year)2014 -2.564e-01 1.740e-01 1.009e+03
## as.numeric(paste(redData$sample_interval)) -7.172e-05 5.697e-05 1.009e+03
## t value Pr(>|t|)
## (Intercept) 1.720 0.085666 .
## redData$AgeGroupyoung -3.880 0.000111 ***
## redData$feed_group2 -0.485 0.627658
## redData$feed_groupNULL -1.157 0.247431
## redData$genetic_groupS -1.383 0.166882
## as.factor(redData$sample_year)2010 -2.013 0.044340 *
## as.factor(redData$sample_year)2011 -1.804 0.071475 .
## as.factor(redData$sample_year)2012 -1.663 0.096544 .
## as.factor(redData$sample_year)2013 -1.467 0.142608
## as.factor(redData$sample_year)2014 -1.474 0.140820
## as.numeric(paste(redData$sample_interval)) -1.259 0.208366
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) rdD$AG rdD$_2 rD$_NU rdD$_S a.(D$_)2010 a.(D$_)2011
## rdDt$AgGrpy -0.106
## rdDt$fd_gr2 -0.069 0.010
## rdDt$f_NULL -0.034 -0.078 0.383
## rdDt$gntc_S -0.060 -0.031 -0.020 0.237
## a.(D$_)2010 -0.983 0.009 0.029 0.013 0.024
## a.(D$_)2011 -0.987 0.033 0.030 0.005 0.022 0.991
## a.(D$_)2012 -0.989 0.044 0.031 0.010 0.024 0.992 0.995
## a.(D$_)2013 -0.989 0.054 0.030 0.006 0.022 0.992 0.995
## a.(D$_)2014 -0.987 0.051 0.027 0.006 0.019 0.991 0.995
## as.n((D$_)) -0.050 0.404 0.044 -0.059 0.056 -0.088 -0.084
## a.(D$_)2012 a.(D$_)2013 a.(D$_)2014
## rdDt$AgGrpy
## rdDt$fd_gr2
## rdDt$f_NULL
## rdDt$gntc_S
## a.(D$_)2010
## a.(D$_)2011
## a.(D$_)2012
## a.(D$_)2013 0.996
## a.(D$_)2014 0.996 0.996
## as.n((D$_)) -0.073 -0.078 -0.090
## convergence code: 0
## boundary (singular) fit: see ?isSingularanova(red_mod_fac)## Type III Analysis of Variance Table with Satterthwaite's method
## Sum Sq Mean Sq NumDF DenDF F value
## redData$AgeGroup 0.44503 0.44503 1 1009 15.0564
## redData$feed_group 0.03965 0.01982 2 1009 0.6707
## redData$genetic_group 0.05656 0.05656 1 1009 1.9135
## as.factor(redData$sample_year) 0.89935 0.17987 5 1009 6.0855
## as.numeric(paste(redData$sample_interval)) 0.04684 0.04684 1 1009 1.5848
## Pr(>F)
## redData$AgeGroup 0.0001111 ***
## redData$feed_group 0.5115971
## redData$genetic_group 0.1668816
## as.factor(redData$sample_year) 1.46e-05 ***
## as.numeric(paste(redData$sample_interval)) 0.2083656
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1AIC(red_mod_fac)## [1] -600.9877anova(red_mod, red_mod_fac)## refitting model(s) with ML (instead of REML)## Data: NULL
## Models:
## red_mod: as.numeric(paste(redData$diff_rltl_res)) ~ redData$age_y + redData$feed_group +
## red_mod: redData$genetic_group + as.factor(redData$sample_year) +
## red_mod: as.numeric(paste(redData$sample_interval)) + (1 | redData$recoded_id)
## red_mod_fac: as.numeric(paste(redData$diff_rltl_res)) ~ redData$AgeGroup +
## red_mod_fac: redData$feed_group + redData$genetic_group + as.factor(redData$sample_year) +
## red_mod_fac: as.numeric(paste(redData$sample_interval)) + (1 | redData$recoded_id)
## npar AIC BIC logLik deviance Chisq Df Pr(>Chisq)
## red_mod 13 -685.40 -621.34 355.70 -711.40
## red_mod_fac 13 -682.27 -618.21 354.14 -708.27 0 0 1Effect of milk productivity on RLTL change
Model of telomere change including milk productivity
#t test with average productivity in the model for animals that have productivity measurements:
PData<-subset(redData, redData$amp_305_lactation != "NULL")
nrow(PData)## [1] 918length(unique(PData$recoded_id))## [1] 253rescaled_prod<-as.numeric(paste(PData$amp_305_lactation))/1000
# milk productivity measures are on a completely differenc scale and are therefore re-scaled by dividing them by 1000
prod_mod<-lmer(as.numeric(paste(PData$diff_rltl_res))~ PData$age_y + PData$feed_group + PData$genetic_group+ as.factor(PData$sample_year) + as.numeric(paste(PData$sample_interval)) + rescaled_prod + (1|PData$recoded_id), na.action=na.exclude)## boundary (singular) fit: see ?isSingularsummary(prod_mod)## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula:
## as.numeric(paste(PData$diff_rltl_res)) ~ PData$age_y + PData$feed_group +
## PData$genetic_group + as.factor(PData$sample_year) + as.numeric(paste(PData$sample_interval)) +
## rescaled_prod + (1 | PData$recoded_id)
##
## REML criterion at convergence: -555.3
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -3.5406 -0.6820 0.0102 0.7086 3.4856
##
## Random effects:
## Groups Name Variance Std.Dev.
## PData$recoded_id (Intercept) 0.00000 0.0000
## Residual 0.02931 0.1712
## Number of obs: 918, groups: PData$recoded_id, 253
##
## Fixed effects:
## Estimate Std. Error df
## (Intercept) 1.954e-01 1.744e-01 9.060e+02
## PData$age_y 2.069e-02 6.722e-03 9.060e+02
## PData$feed_group2 -5.600e-03 1.305e-02 9.060e+02
## PData$feed_groupNULL 9.214e-03 2.429e-02 9.060e+02
## PData$genetic_groupS -1.131e-02 1.222e-02 9.060e+02
## as.factor(PData$sample_year)2010 -3.381e-01 1.739e-01 9.060e+02
## as.factor(PData$sample_year)2011 -3.005e-01 1.737e-01 9.060e+02
## as.factor(PData$sample_year)2012 -2.783e-01 1.736e-01 9.060e+02
## as.factor(PData$sample_year)2013 -2.522e-01 1.745e-01 9.060e+02
## as.factor(PData$sample_year)2014 -2.746e-01 1.757e-01 9.060e+02
## as.numeric(paste(PData$sample_interval)) -2.573e-05 5.582e-05 9.060e+02
## rescaled_prod 2.735e-03 2.840e-03 9.060e+02
## t value Pr(>|t|)
## (Intercept) 1.121 0.26269
## PData$age_y 3.078 0.00215 **
## PData$feed_group2 -0.429 0.66799
## PData$feed_groupNULL 0.379 0.70449
## PData$genetic_groupS -0.926 0.35490
## as.factor(PData$sample_year)2010 -1.945 0.05214 .
## as.factor(PData$sample_year)2011 -1.730 0.08395 .
## as.factor(PData$sample_year)2012 -1.603 0.10931
## as.factor(PData$sample_year)2013 -1.446 0.14862
## as.factor(PData$sample_year)2014 -1.563 0.11845
## as.numeric(paste(PData$sample_interval)) -0.461 0.64498
## rescaled_prod 0.963 0.33571
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) PDt$g_ PDt$_2 PD$_NU PDt$_S a.(PD$_)2010 a.(PD$_)2011
## PData$age_y 0.033
## PDt$fd_grp2 -0.124 -0.163
## PDt$fd_NULL -0.127 -0.143 0.452
## PDt$gntc_gS -0.016 0.102 -0.130 -0.006
## a.(PD$_)2010 -0.980 -0.049 0.050 0.044 0.010
## a.(PD$_)2011 -0.981 -0.069 0.054 0.044 0.006 0.991
## a.(PD$_)2012 -0.982 -0.087 0.057 0.051 0.005 0.991 0.994
## a.(PD$_)2013 -0.979 -0.120 0.062 0.056 0.000 0.990 0.994
## a.(PD$_)2014 -0.972 -0.161 0.065 0.060 -0.004 0.986 0.990
## as.((PD$_)) -0.001 0.008 0.033 -0.064 0.059 -0.112 -0.115
## rescald_prd -0.164 -0.178 0.376 0.599 -0.283 0.055 0.057
## a.(PD$_)2012 a.(PD$_)2013 a.(PD$_)2014 a.((PD
## PData$age_y
## PDt$fd_grp2
## PDt$fd_NULL
## PDt$gntc_gS
## a.(PD$_)2010
## a.(PD$_)2011
## a.(PD$_)2012
## a.(PD$_)2013 0.995
## a.(PD$_)2014 0.993 0.995
## as.((PD$_)) -0.110 -0.121 -0.131
## rescald_prd 0.060 0.068 0.071 -0.036
## convergence code: 0
## boundary (singular) fit: see ?isSingularanova(prod_mod)## Type III Analysis of Variance Table with Satterthwaite's method
## Sum Sq Mean Sq NumDF DenDF F value
## PData$age_y 0.27757 0.277574 1 906 9.4712
## PData$feed_group 0.01751 0.008753 2 906 0.2987
## PData$genetic_group 0.02511 0.025108 1 906 0.8567
## as.factor(PData$sample_year) 0.49716 0.099433 5 906 3.3928
## as.numeric(paste(PData$sample_interval)) 0.00623 0.006226 1 906 0.2124
## rescaled_prod 0.02719 0.027189 1 906 0.9277
## Pr(>F)
## PData$age_y 0.002150 **
## PData$feed_group 0.741877
## PData$genetic_group 0.354903
## as.factor(PData$sample_year) 0.004824 **
## as.numeric(paste(PData$sample_interval)) 0.644976
## rescaled_prod 0.335712
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Analysis of change between first two measurements within first year and productive lifespan
Histogram of change in telomere length within the first year of life
ggplot(init_change_dat, aes(x=init_change_dat$ResDiff)) +geom_histogram(colour= "black", fill = mycoloursP[6])+ xlab("Initial RLTL change") + theme_classic(20)+
geom_vline(xintercept= mean(as.numeric(paste(init_change_dat$ResDiff))), colour="red") +
annotate("text", x = -0.7, y = 28, label = "(*)", fontface =2, size= 10) +
stat_function(
fun = function(x, mean, sd, n){
n/23 * dnorm(x = x, mean = mean, sd = sd)
},
args = with(init_change_dat, c(mean = mean(as.numeric(paste(init_change_dat$ResDiff))),
sd = sd(as.numeric(paste(init_change_dat$ResDiff))),
n= length(as.numeric(paste(init_change_dat$ResDiff))))))## Warning: Use of `init_change_dat$ResDiff` is discouraged. Use `ResDiff` instead.
## Warning: Use of `init_change_dat$ResDiff` is discouraged. Use `ResDiff` instead.## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
Statistical normality tests
shapiro.test(init_change_dat$ResDiff)##
## Shapiro-Wilk normality test
##
## data: init_change_dat$ResDiff
## W = 0.99401, p-value = 0.3057ks.test(init_change_dat$ResDiff, pnorm)##
## One-sample Kolmogorov-Smirnov test
##
## data: init_change_dat$ResDiff
## D = 0.39675, p-value < 2.2e-16
## alternative hypothesis: two-sidedqqnorm(init_change_dat$ResDiff)
Cox-proportional hazard model of productive lifespan
coxFit<-coxph(Surv(as.numeric(paste(init_change_dat$TimeMeasure)), init_change_dat$Event == 1)~init_change_dat$ResDiff)
summary(coxFit)## Call:
## coxph(formula = Surv(as.numeric(paste(init_change_dat$TimeMeasure)),
## init_change_dat$Event == 1) ~ init_change_dat$ResDiff)
##
## n= 291, number of events= 227
##
## coef exp(coef) se(coef) z Pr(>|z|)
## init_change_dat$ResDiff -1.1406 0.3196 0.3914 -2.914 0.00356 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## exp(coef) exp(-coef) lower .95 upper .95
## init_change_dat$ResDiff 0.3196 3.129 0.1484 0.6883
##
## Concordance= 0.554 (se = 0.021 )
## Likelihood ratio test= 8.32 on 1 df, p=0.004
## Wald test = 8.49 on 1 df, p=0.004
## Score (logrank) test = 8.48 on 1 df, p=0.004change telomere change measure to discrete scale to allow for plotting
ggplot(init_change_dat, aes(x=as.factor(init_change_dat$ResDiff_3G), y=init_change_dat$ResDiff)) + geom_boxplot() +xlab("RLTL change tertiles") + ylab("RLTL change")+ theme_classic(20)+ geom_jitter()## Warning: Use of `init_change_dat$ResDiff_3G` is discouraged. Use `ResDiff_3G`
## instead.## Warning: Use of `init_change_dat$ResDiff` is discouraged. Use `ResDiff` instead.## Warning: Use of `init_change_dat$ResDiff_3G` is discouraged. Use `ResDiff_3G`
## instead.## Warning: Use of `init_change_dat$ResDiff` is discouraged. Use `ResDiff` instead.
Cox proportional hazard model on discrete scale change measurements
coxFit<-coxph(Surv(init_change_dat$TimeMeasure, init_change_dat$Event == 1)~init_change_dat$ResDiff_3G)
summary(coxFit)## Call:
## coxph(formula = Surv(init_change_dat$TimeMeasure, init_change_dat$Event ==
## 1) ~ init_change_dat$ResDiff_3G)
##
## n= 291, number of events= 227
##
## coef exp(coef) se(coef) z Pr(>|z|)
## init_change_dat$ResDiff_3G -0.2259 0.7978 0.0820 -2.755 0.00586 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## exp(coef) exp(-coef) lower .95 upper .95
## init_change_dat$ResDiff_3G 0.7978 1.253 0.6793 0.9369
##
## Concordance= 0.55 (se = 0.02 )
## Likelihood ratio test= 7.61 on 1 df, p=0.006
## Wald test = 7.59 on 1 df, p=0.006
## Score (logrank) test = 7.65 on 1 df, p=0.006Plot discrete measures
#g1<- ggsurvplot(coxFit, data=init_change_dat, pval = 0.006, ggtheme = theme_minimal(), risk.table.y.text.col = T, risk.table.y.text = FALSE, legend.labs=c("Tertile 1 (most RLTL attrition)", "Tertile 2", "Tertile 3 (least RLTL attrition"), surv.median.line = "hv")+ guides(colour = guide_legend(nrow = 3))+ theme_survminer(40)
#g1
#ggsurvplot(coxFit, data=init_change_dat, pval = 0.006, ggtheme = theme_minimal(), risk.table.y.text.col = T, risk.table.y.text = FALSE, legend.labs=c("Tertile 1 (most RLTL attrition)", "Tertile 2", "Tertile 3 (least RLTL attrition"), surv.median.line = "hv")+theme_survminer(40) + xlab("Time in days")
#ggsave(file = "/Users/luise/Documents/Work/CattleTelomereLength/AgingCell/Images/InitialChange.pdf", print(g1))Individual longitudinal RLTL data of animals with more than seven measurements
spaghettiData<-read.csv("/Users/luise/Documents/Work/CattleTelomereLength/LongTelomereChangePaper/git/TelomereChangeInDairyCattle/RC_SpaghettiData.csv")
ggplot(data = spaghettiData, aes(x = age_y, y = residual_rltl_t, group = recoded_id))+geom_point()+ geom_line()+facet_grid(.~recoded_id)+
xlab("Age in years") + ylab("RLTL residuals")
ggplot(data = data, aes(x = sample_year, y = residual_rltl_t, group = recoded_id, colour= as.factor(data$recoded_id)))+geom_point()+ geom_line()+
xlab("Age in years") + ylab("RLTL residuals") +
theme(legend.position = "none")## Warning: Use of `data$recoded_id` is discouraged. Use `recoded_id` instead.
## Warning: Use of `data$recoded_id` is discouraged. Use `recoded_id` instead.
ggplot(data = data, aes(x = age_y , y = residual_rltl_t, group = recoded_id)) + geom_point(color= "white")+ geom_line(color= "white") + xlab("Age in years") + ylab("RLTL residuals") + theme(legend.position = "none") + geom_smooth(method = "lm", se = FALSE) + theme_classic()## `geom_smooth()` using formula 'y ~ x'
# # average_resid rltl_change abs_rltl_change MeanChangeRate MeanChangeRateAbs
deadD<-subset(data, data$herd_life != "NULL")
redDeadD<-subset(redData, redData$herd_life != "NULL")
uDeadData<-redDeadD[!duplicated(redDeadD$recoded_id),]
#1)
x<-15
AverRes<-uDeadData$average_resid
chF<-uDeadData$rltl_change
df1<-data.frame(AverRes, chF)
cor.test(AverRes, chF)##
## Pearson's product-moment correlation
##
## data: AverRes and chF
## t = -3.0223, df = 239, p-value = 0.002782
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
## -0.31069811 -0.06712565
## sample estimates:
## cor
## -0.1918646CorP1 <- ggplot(df1 , aes(x= AverRes, y=chF)) + geom_point(size= 2) +
labs(x= "Mean RLTL", y= "Mean RLTL change") + stat_smooth(method=lm, formula = y~x)+
theme_gray(x)
CorP1
#2)
AchF<-uDeadData$abs_rltl_change
df2<-data.frame(AverRes, AchF)
cor.test(AverRes, AchF)##
## Pearson's product-moment correlation
##
## data: AverRes and AchF
## t = 4.9383, df = 239, p-value = 1.48e-06
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
## 0.185034 0.414706
## sample estimates:
## cor
## 0.3042856CorP2 <- ggplot(df2 , aes(x= AverRes, y=AchF)) + geom_point(size= 2) +
labs(x= "Mean RLTL", y= "Mean absolute RLTL change") + stat_smooth(method=lm, formula = y~x)+
theme_gray(x)
CorP2
#5)
df5<-data.frame(chF, AchF)
cor.test(chF, AchF)##
## Pearson's product-moment correlation
##
## data: chF and AchF
## t = -9.7941, df = 239, p-value < 2.2e-16
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
## -0.6196323 -0.4384552
## sample estimates:
## cor
## -0.53517CorP5 <- ggplot(df5 , aes(x= chF, y=AchF)) + geom_point(size= 2) +
labs(x= "Mean RLTL change", y= "Mean absolute RLTL change") +
theme_gray(x)+ stat_smooth(method=lm, formula = y~x)+
geom_vline(xintercept=0, linetype="dotted") +
geom_abline(intercept=0, slope=-1, colour="red")+
geom_abline(intercept=0, slope=1, colour= "blue")
CorP5
#grid.arrange(CorP1, CorP2, CorP3, CorP4, CorP5, CorP6, CorP7, CorP8, CorP9, CorP10, ncol=3 ) other plots can be found in 20170914_ChangeAnalysis
grid.arrange(CorP1, CorP2, CorP5, ncol=3 )
Correlation between subsquent measurements
cor.test(redData$residual_rltl_t,as.numeric(paste(redData$residual_rltl_tminus1)))##
## Pearson's product-moment correlation
##
## data: redData$residual_rltl_t and as.numeric(paste(redData$residual_rltl_tminus1))
## t = 13.058, df = 1018, p-value < 2.2e-16
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
## 0.3249385 0.4301473
## sample estimates:
## cor
## 0.3787659SubsPlot<-ggplot(redData, aes(x=as.numeric(paste(residual_rltl_tminus1)), y=redData$residual_rltl_t))+geom_point()+xlab("RLTL at time t-1")+ylab("RLTL at time t")+theme_gray(20)+geom_abline(intercept=0, slope=1, colour="red")+ stat_smooth(method=lm, formula=y~x)
SubsPlot## Warning: Use of `redData$residual_rltl_t` is discouraged. Use `residual_rltl_t`
## instead.
## Warning: Use of `redData$residual_rltl_t` is discouraged. Use `residual_rltl_t`
## instead.
#Plot correlation between present and past change
redRedData<-subset(data, data$change_at_t_1!= "NULL")
nrow(redRedData)## [1] 715cor.test(as.numeric(paste(redRedData$diff_rltl_res)), as.numeric(paste(redRedData$change_at_t_1)))##
## Pearson's product-moment correlation
##
## data: as.numeric(paste(redRedData$diff_rltl_res)) and as.numeric(paste(redRedData$change_at_t_1))
## t = -12.838, df = 713, p-value < 2.2e-16
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
## -0.4910345 -0.3718059
## sample estimates:
## cor
## -0.4333142SubsPlot2<-ggplot(redRedData, aes(x=as.numeric(paste(redRedData$change_at_t_1)), y=as.numeric(paste(redRedData$diff_rltl_res))))+geom_point()+xlab("RLTL change at time t-1")+ylab("RLTL change at time t")+theme_gray(20)+geom_abline(intercept=0, slope=-1, colour="red")+ stat_smooth(method=lm, formula=y~x)
SubsPlot2## Warning: Use of `redRedData$change_at_t_1` is discouraged. Use `change_at_t_1`
## instead.## Warning: Use of `redRedData$diff_rltl_res` is discouraged. Use `diff_rltl_res`
## instead.## Warning: Use of `redRedData$change_at_t_1` is discouraged. Use `change_at_t_1`
## instead.## Warning: Use of `redRedData$diff_rltl_res` is discouraged. Use `diff_rltl_res`
## instead.
grid.arrange(SubsPlot, SubsPlot2, ncol=2)## Warning: Use of `redData$residual_rltl_t` is discouraged. Use `residual_rltl_t`
## instead.## Warning: Use of `redData$residual_rltl_t` is discouraged. Use `residual_rltl_t`
## instead.## Warning: Use of `redRedData$change_at_t_1` is discouraged. Use `change_at_t_1`
## instead.## Warning: Use of `redRedData$diff_rltl_res` is discouraged. Use `diff_rltl_res`
## instead.## Warning: Use of `redRedData$change_at_t_1` is discouraged. Use `change_at_t_1`
## instead.## Warning: Use of `redRedData$diff_rltl_res` is discouraged. Use `diff_rltl_res`
## instead.
#resRTLatTmin1<-as.numeric(paste(redData$residual_rltl_tminus1))Productive lifespan
variable<- as.numeric(paste(uDeadData$herd_life))
#sqrt might improve normality test outcome slightly...
meanHL<-mean(as.numeric(paste(uDeadData$herd_life)))
prod_lsp_plot <- ggplot(uDeadData, aes(x = as.numeric(paste(uDeadData$herd_life)))) +
geom_histogram(colour= "black", fill = mycoloursP[4]) + xlab("Productive lifespan in days") +
geom_vline(xintercept= mean(as.numeric(paste(uDeadData$herd_life))), colour="red")+
annotate("text", x = 0, y = 17.5, label = "(D)", fontface =2, size= 10) + theme_classic(20) +
stat_function(
fun = function(x, mean, sd, n){
n*120 * dnorm(x = x, mean = mean, sd = sd)
},
args = with(uDeadData, c(mean = mean(as.numeric(paste(uDeadData$herd_life))),
sd = sd(as.numeric(paste(uDeadData$herd_life))),
n= length(as.numeric(paste(uDeadData$herd_life))))))
prod_lsp_plot## Warning: Use of `uDeadData$herd_life` is discouraged. Use `herd_life` instead.
## Warning: Use of `uDeadData$herd_life` is discouraged. Use `herd_life` instead.## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
# statistical tests for normal distribution of RTL measurements
shapiro.test(variable)##
## Shapiro-Wilk normality test
##
## data: variable
## W = 0.98408, p-value = 0.008487ks.test(variable, pnorm)## Warning in ks.test(variable, pnorm): ties should not be present for the
## Kolmogorov-Smirnov test##
## One-sample Kolmogorov-Smirnov test
##
## data: variable
## D = 1, p-value < 2.2e-16
## alternative hypothesis: two-sidedqqnorm(variable)
Assemble plots for supplemenrary file
grid.arrange(hist1, hist2, hist3, prod_lsp_plot, ncol=2)## Warning: Use of `uData$SamplesPerAnimal` is discouraged. Use `SamplesPerAnimal`
## instead.## Warning: Use of `data$residual_rltl_t` is discouraged. Use `residual_rltl_t`
## instead.
## Warning: Use of `data$residual_rltl_t` is discouraged. Use `residual_rltl_t`
## instead.## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.## Warning: Use of `redData$diff_rltl_res` is discouraged. Use `diff_rltl_res`
## instead.## Warning: Use of `redData$diff_rltl_res` is discouraged. Use `diff_rltl_res`
## instead.## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.## Warning: Use of `uDeadData$herd_life` is discouraged. Use `herd_life` instead.## Warning: Use of `uDeadData$herd_life` is discouraged. Use `herd_life` instead.## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
Histograms of (A) Number of RLTL measurements per animal, (B) RLTL residuals after correcting for qPCR plate and row, (C) measurement to measurement change in RLTL residuals (after correcting for qPCR plate and row) and (D) productive lifespan in days; Red lines indicate mean.
Backwards elimination of all non-significant effects
#Age in years as factor with reference year being 0 (Age at t-1)
mod1f<-lmer(as.numeric(paste(redData$diff_rltl_res))~ as.factor(redData$age_y) + (1|redData$recoded_id), na.action=na.exclude)## boundary (singular) fit: see ?isSingularsummary(mod1f)## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: as.numeric(paste(redData$diff_rltl_res)) ~ as.factor(redData$age_y) +
## (1 | redData$recoded_id)
##
## REML criterion at convergence: -691.2
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -3.6127 -0.6744 0.0379 0.6816 3.5134
##
## Random effects:
## Groups Name Variance Std.Dev.
## redData$recoded_id (Intercept) 0.00000 0.0000
## Residual 0.02871 0.1694
## Number of obs: 1020, groups: redData$recoded_id, 305
##
## Fixed effects:
## Estimate Std. Error df t value Pr(>|t|)
## (Intercept) -0.10057 0.06404 1013.00000 -1.570 0.1167
## as.factor(redData$age_y)1 -0.01458 0.06482 1013.00000 -0.225 0.8221
## as.factor(redData$age_y)2 0.09171 0.06484 1013.00000 1.414 0.1575
## as.factor(redData$age_y)3 0.11267 0.06503 1013.00000 1.733 0.0835 .
## as.factor(redData$age_y)4 0.09686 0.06558 1013.00000 1.477 0.1400
## as.factor(redData$age_y)5 0.12766 0.06800 1013.00000 1.877 0.0608 .
## as.factor(redData$age_y)6 0.07425 0.07441 1013.00000 0.998 0.3186
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) a.(D$_)1 a.(D$_)2 a.(D$_)3 a.(D$_)4 a.(D$_)5
## as.fc(D$_)1 -0.988
## as.fc(D$_)2 -0.988 0.976
## as.fc(D$_)3 -0.985 0.973 0.973
## as.fc(D$_)4 -0.977 0.965 0.965 0.962
## as.fc(D$_)5 -0.942 0.931 0.930 0.928 0.920
## as.fc(D$_)6 -0.861 0.850 0.850 0.848 0.840 0.811
## convergence code: 0
## boundary (singular) fit: see ?isSingularanova(mod1f)## Type III Analysis of Variance Table with Satterthwaite's method
## Sum Sq Mean Sq NumDF DenDF F value Pr(>F)
## as.factor(redData$age_y) 2.9149 0.48582 6 1013 16.921 < 2.2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1################
################
#Test
#mod1f<-lmer(abs(as.numeric(paste(redData$diff_rltl_res)))~ 1 + (1|redData$recoded_id), na.action=na.exclude)
#summary(mod1f)
#anova(mod1f)
0.0006458/(0.0006458+0.0111165)## [1] 0.05490423##############
###############
redData$R_age_y<-as.factor(redData$age_y)
#lm without random effect to do likelihood ratio test
modLM<-glm(as.numeric(paste(redData$diff_rltl_res))~ as.factor(redData$age_y))
anova(mod1f, modLM)## refitting model(s) with ML (instead of REML)## Data: NULL
## Models:
## modLM: as.numeric(paste(redData$diff_rltl_res)) ~ as.factor(redData$age_y)
## mod1f: as.numeric(paste(redData$diff_rltl_res)) ~ as.factor(redData$age_y) +
## mod1f: (1 | redData$recoded_id)
## npar AIC BIC logLik deviance Chisq Df Pr(>Chisq)
## modLM 8 -717.85 -678.43 366.93 -733.85
## mod1f 9 -715.85 -671.50 366.93 -733.85 0 1 1logLik(mod1f)## 'log Lik.' 345.6007 (df=9)logLik(modLM)## 'log Lik.' 366.9254 (df=8)#Age in years as factor with reference year being 1 (Age at t-1)
contrasts(redData$R_age_y) <- contr.treatment(levels(redData$R_age_y), base=which(levels(redData$R_age_y) == '1'))
mod1f<-lmer(as.numeric(paste(redData$diff_rltl_res))~ as.factor(redData$R_age_y)+ as.factor(redData$sample_year) + as.numeric(paste(redData$sample_interval)) + (1|redData$recoded_id), na.action=na.exclude)## boundary (singular) fit: see ?isSingularsummary(mod1f)## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula:
## as.numeric(paste(redData$diff_rltl_res)) ~ as.factor(redData$R_age_y) +
## as.factor(redData$sample_year) + as.numeric(paste(redData$sample_interval)) +
## (1 | redData$recoded_id)
##
## REML criterion at convergence: -653.9
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -3.5347 -0.6741 0.0278 0.6901 3.5959
##
## Random effects:
## Groups Name Variance Std.Dev.
## redData$recoded_id (Intercept) 0.00000 0.0000
## Residual 0.02868 0.1693
## Number of obs: 1020, groups: redData$recoded_id, 305
##
## Fixed effects:
## Estimate Std. Error df
## (Intercept) 2.363e-01 1.842e-01 1.007e+03
## as.factor(redData$R_age_y)0 -3.155e-02 7.258e-02 1.007e+03
## as.factor(redData$R_age_y)2 9.681e-02 1.561e-02 1.007e+03
## as.factor(redData$R_age_y)3 1.131e-01 1.827e-02 1.007e+03
## as.factor(redData$R_age_y)4 9.575e-02 2.220e-02 1.007e+03
## as.factor(redData$R_age_y)5 1.322e-01 3.037e-02 1.007e+03
## as.factor(redData$R_age_y)6 7.947e-02 4.352e-02 1.007e+03
## as.factor(redData$sample_year)2010 -3.634e-01 1.831e-01 1.007e+03
## as.factor(redData$sample_year)2011 -3.539e-01 1.839e-01 1.007e+03
## as.factor(redData$sample_year)2012 -3.420e-01 1.838e-01 1.007e+03
## as.factor(redData$sample_year)2013 -3.257e-01 1.838e-01 1.007e+03
## as.factor(redData$sample_year)2014 -3.409e-01 1.839e-01 1.007e+03
## as.numeric(paste(redData$sample_interval)) -3.703e-06 5.319e-05 1.007e+03
## t value Pr(>|t|)
## (Intercept) 1.282 0.2000
## as.factor(redData$R_age_y)0 -0.435 0.6639
## as.factor(redData$R_age_y)2 6.201 8.19e-10 ***
## as.factor(redData$R_age_y)3 6.192 8.64e-10 ***
## as.factor(redData$R_age_y)4 4.314 1.76e-05 ***
## as.factor(redData$R_age_y)5 4.352 1.49e-05 ***
## as.factor(redData$R_age_y)6 1.826 0.0681 .
## as.factor(redData$sample_year)2010 -1.984 0.0475 *
## as.factor(redData$sample_year)2011 -1.924 0.0546 .
## as.factor(redData$sample_year)2012 -1.861 0.0631 .
## as.factor(redData$sample_year)2013 -1.772 0.0768 .
## as.factor(redData$sample_year)2014 -1.854 0.0641 .
## as.numeric(paste(redData$sample_interval)) -0.070 0.9445
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1##
## Correlation matrix not shown by default, as p = 13 > 12.
## Use print(x, correlation=TRUE) or
## vcov(x) if you need it## convergence code: 0
## boundary (singular) fit: see ?isSingularanova(mod1f)## Type III Analysis of Variance Table with Satterthwaite's method
## Sum Sq Mean Sq NumDF DenDF F value
## as.factor(redData$R_age_y) 1.49717 0.249528 6 1007 8.7016
## as.factor(redData$sample_year) 0.20720 0.041441 5 1007 1.4451
## as.numeric(paste(redData$sample_interval)) 0.00014 0.000139 1 1007 0.0048
## Pr(>F)
## as.factor(redData$R_age_y) 2.959e-09 ***
## as.factor(redData$sample_year) 0.2054
## as.numeric(paste(redData$sample_interval)) 0.9445
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1exclude2009<-subset(redData, redData$sample_year != "2009")
mod1f<-lmer(as.numeric(paste(exclude2009$diff_rltl_res))~ as.factor(exclude2009$R_age_y)+ as.factor(exclude2009$sample_year) + as.numeric(paste(exclude2009$sample_interval)) + (1|exclude2009$recoded_id), na.action=na.exclude)## boundary (singular) fit: see ?isSingularsummary(mod1f)## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula:
## as.numeric(paste(exclude2009$diff_rltl_res)) ~ as.factor(exclude2009$R_age_y) +
## as.factor(exclude2009$sample_year) + as.numeric(paste(exclude2009$sample_interval)) +
## (1 | exclude2009$recoded_id)
##
## REML criterion at convergence: -653.9
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -3.5347 -0.6744 0.0286 0.6916 3.5959
##
## Random effects:
## Groups Name Variance Std.Dev.
## exclude2009$recoded_id (Intercept) 0.00000 0.0000
## Residual 0.02868 0.1693
## Number of obs: 1019, groups: exclude2009$recoded_id, 305
##
## Fixed effects:
## Estimate Std. Error df
## (Intercept) -1.272e-01 2.672e-02 1.007e+03
## as.factor(exclude2009$R_age_y)0 -3.155e-02 7.258e-02 1.007e+03
## as.factor(exclude2009$R_age_y)2 9.681e-02 1.561e-02 1.007e+03
## as.factor(exclude2009$R_age_y)3 1.131e-01 1.827e-02 1.007e+03
## as.factor(exclude2009$R_age_y)4 9.575e-02 2.220e-02 1.007e+03
## as.factor(exclude2009$R_age_y)5 1.322e-01 3.037e-02 1.007e+03
## as.factor(exclude2009$R_age_y)6 7.947e-02 4.352e-02 1.007e+03
## as.factor(exclude2009$sample_year)2011 9.552e-03 2.315e-02 1.007e+03
## as.factor(exclude2009$sample_year)2012 2.146e-02 2.309e-02 1.007e+03
## as.factor(exclude2009$sample_year)2013 3.772e-02 2.459e-02 1.007e+03
## as.factor(exclude2009$sample_year)2014 2.254e-02 2.724e-02 1.007e+03
## as.numeric(paste(exclude2009$sample_interval)) -3.703e-06 5.319e-05 1.007e+03
## t value Pr(>|t|)
## (Intercept) -4.759 2.23e-06 ***
## as.factor(exclude2009$R_age_y)0 -0.435 0.6639
## as.factor(exclude2009$R_age_y)2 6.201 8.19e-10 ***
## as.factor(exclude2009$R_age_y)3 6.192 8.64e-10 ***
## as.factor(exclude2009$R_age_y)4 4.314 1.76e-05 ***
## as.factor(exclude2009$R_age_y)5 4.352 1.49e-05 ***
## as.factor(exclude2009$R_age_y)6 1.826 0.0681 .
## as.factor(exclude2009$sample_year)2011 0.413 0.6799
## as.factor(exclude2009$sample_year)2012 0.930 0.3528
## as.factor(exclude2009$sample_year)2013 1.534 0.1253
## as.factor(exclude2009$sample_year)2014 0.827 0.4083
## as.numeric(paste(exclude2009$sample_interval)) -0.070 0.9445
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) a.(2009$R__)0 a.(2009$R__)2 a.(2009$R__)3 a.(2009$R__)4
## a.(2009$R__)0 -0.287
## a.(2009$R__)2 -0.031 0.057
## a.(2009$R__)3 -0.007 0.064 0.560
## a.(2009$R__)4 0.019 0.060 0.495 0.585
## a.(2009$R__)5 -0.015 0.066 0.362 0.456 0.483
## a.(2009$R__)6 -0.107 0.078 0.251 0.318 0.335
## a.(2009$_)2011 -0.520 0.113 -0.267 -0.162 -0.133
## a.(2009$_)2012 -0.566 0.117 -0.329 -0.336 -0.236
## a.(2009$_)2013 -0.486 0.081 -0.405 -0.453 -0.433
## a.(2009$_)2014 -0.362 0.013 -0.362 -0.484 -0.516
## a.((2009$_) -0.726 0.232 -0.014 -0.028 -0.062
## a.(2009$R__)5 a.(2009$R__)6 a.(2009$_)2011 a.(2009$_)2012
## a.(2009$R__)0
## a.(2009$R__)2
## a.(2009$R__)3
## a.(2009$R__)4
## a.(2009$R__)5
## a.(2009$R__)6 0.313
## a.(2009$_)2011 -0.097 -0.070
## a.(2009$_)2012 -0.173 -0.115 0.719
## a.(2009$_)2013 -0.281 -0.192 0.702 0.775
## a.(2009$_)2014 -0.480 -0.362 0.632 0.708
## a.((2009$_) -0.008 0.127 -0.017 0.049
## a.(2009$_)2013 a.(2009$_)2014
## a.(2009$R__)0
## a.(2009$R__)2
## a.(2009$R__)3
## a.(2009$R__)4
## a.(2009$R__)5
## a.(2009$R__)6
## a.(2009$_)2011
## a.(2009$_)2012
## a.(2009$_)2013
## a.(2009$_)2014 0.790
## a.((2009$_) -0.010 -0.108
## convergence code: 0
## boundary (singular) fit: see ?isSingularanova(mod1f)## Type III Analysis of Variance Table with Satterthwaite's method
## Sum Sq Mean Sq NumDF DenDF
## as.factor(exclude2009$R_age_y) 1.49717 0.249528 6 1007
## as.factor(exclude2009$sample_year) 0.09885 0.024712 4 1007
## as.numeric(paste(exclude2009$sample_interval)) 0.00014 0.000139 1 1007
## F value Pr(>F)
## as.factor(exclude2009$R_age_y) 8.7016 2.959e-09 ***
## as.factor(exclude2009$sample_year) 0.8618 0.4864
## as.numeric(paste(exclude2009$sample_interval)) 0.0048 0.9445
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1df<-as.data.frame(round(coef(summary(mod1f)),3))
#write.csv(df, "/Users/luiseseeker/Desktop/PNAS/S2_Table_AgeFactor.csv")
#lsmeans(mod1f,pairwise~as.factor(redData$R_age_y), adjust="Bonferroni")
#Age in years as factor with reference year being 2 (Age at t-1)
contrasts(redData$R_age_y) <- contr.treatment(levels(redData$R_age_y), base=which(levels(redData$R_age_y) == '2'))
mod1f<-lmer(as.numeric(paste(redData$diff_rltl_res))~ as.factor(redData$R_age_y) + (1|redData$recoded_id), na.action=na.exclude)## boundary (singular) fit: see ?isSingularsummary(mod1f)## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula:
## as.numeric(paste(redData$diff_rltl_res)) ~ as.factor(redData$R_age_y) +
## (1 | redData$recoded_id)
##
## REML criterion at convergence: -691.2
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -3.6127 -0.6744 0.0379 0.6816 3.5134
##
## Random effects:
## Groups Name Variance Std.Dev.
## redData$recoded_id (Intercept) 0.00000 0.0000
## Residual 0.02871 0.1694
## Number of obs: 1020, groups: redData$recoded_id, 305
##
## Fixed effects:
## Estimate Std. Error df t value Pr(>|t|)
## (Intercept) -8.857e-03 1.011e-02 1.013e+03 -0.876 0.381
## as.factor(redData$R_age_y)0 -9.171e-02 6.484e-02 1.013e+03 -1.414 0.158
## as.factor(redData$R_age_y)1 -1.063e-01 1.422e-02 1.013e+03 -7.474 1.67e-13
## as.factor(redData$R_age_y)3 2.096e-02 1.514e-02 1.013e+03 1.385 0.166
## as.factor(redData$R_age_y)4 5.150e-03 1.737e-02 1.013e+03 0.297 0.767
## as.factor(redData$R_age_y)5 3.595e-02 2.498e-02 1.013e+03 1.439 0.150
## as.factor(redData$R_age_y)6 -1.746e-02 3.921e-02 1.013e+03 -0.445 0.656
##
## (Intercept)
## as.factor(redData$R_age_y)0
## as.factor(redData$R_age_y)1 ***
## as.factor(redData$R_age_y)3
## as.factor(redData$R_age_y)4
## as.factor(redData$R_age_y)5
## as.factor(redData$R_age_y)6
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) a.(D$R__)0 a.(D$R__)1 a.(D$R__)3 a.(D$R__)4 a.(D$R__)5
## as.(D$R__)0 -0.156
## as.(D$R__)1 -0.711 0.111
## as.(D$R__)3 -0.668 0.104 0.475
## as.(D$R__)4 -0.582 0.091 0.414 0.389
## as.(D$R__)5 -0.405 0.063 0.288 0.270 0.236
## as.(D$R__)6 -0.258 0.040 0.183 0.172 0.150 0.104
## convergence code: 0
## boundary (singular) fit: see ?isSingularanova(mod1f)## Type III Analysis of Variance Table with Satterthwaite's method
## Sum Sq Mean Sq NumDF DenDF F value Pr(>F)
## as.factor(redData$R_age_y) 2.9149 0.48582 6 1013 16.921 < 2.2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1#Age in years as factor with reference year being 3 (Age at t-1)
contrasts(redData$R_age_y) <- contr.treatment(levels(redData$R_age_y), base=which(levels(redData$R_age_y) == '3'))
mod1f<-lmer(as.numeric(paste(redData$diff_rltl_res))~ as.factor(redData$R_age_y) + (1|redData$recoded_id), na.action=na.exclude)## boundary (singular) fit: see ?isSingularsummary(mod1f)## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula:
## as.numeric(paste(redData$diff_rltl_res)) ~ as.factor(redData$R_age_y) +
## (1 | redData$recoded_id)
##
## REML criterion at convergence: -691.2
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -3.6127 -0.6744 0.0379 0.6816 3.5134
##
## Random effects:
## Groups Name Variance Std.Dev.
## redData$recoded_id (Intercept) 0.00000 0.0000
## Residual 0.02871 0.1694
## Number of obs: 1020, groups: redData$recoded_id, 305
##
## Fixed effects:
## Estimate Std. Error df t value Pr(>|t|)
## (Intercept) 0.01211 0.01127 1013.00000 1.074 0.2830
## as.factor(redData$R_age_y)0 -0.11267 0.06503 1013.00000 -1.733 0.0835
## as.factor(redData$R_age_y)1 -0.12725 0.01507 1013.00000 -8.444 <2e-16
## as.factor(redData$R_age_y)2 -0.02096 0.01514 1013.00000 -1.385 0.1665
## as.factor(redData$R_age_y)4 -0.01581 0.01807 1013.00000 -0.875 0.3816
## as.factor(redData$R_age_y)5 0.01499 0.02548 1013.00000 0.588 0.5565
## as.factor(redData$R_age_y)6 -0.03843 0.03953 1013.00000 -0.972 0.3312
##
## (Intercept)
## as.factor(redData$R_age_y)0 .
## as.factor(redData$R_age_y)1 ***
## as.factor(redData$R_age_y)2
## as.factor(redData$R_age_y)4
## as.factor(redData$R_age_y)5
## as.factor(redData$R_age_y)6
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) a.(D$R__)0 a.(D$R__)1 a.(D$R__)2 a.(D$R__)4 a.(D$R__)5
## as.(D$R__)0 -0.173
## as.(D$R__)1 -0.748 0.130
## as.(D$R__)2 -0.744 0.129 0.557
## as.(D$R__)4 -0.624 0.108 0.467 0.464
## as.(D$R__)5 -0.442 0.077 0.331 0.329 0.276
## as.(D$R__)6 -0.285 0.049 0.213 0.212 0.178 0.126
## convergence code: 0
## boundary (singular) fit: see ?isSingularanova(mod1f)## Type III Analysis of Variance Table with Satterthwaite's method
## Sum Sq Mean Sq NumDF DenDF F value Pr(>F)
## as.factor(redData$R_age_y) 2.9149 0.48582 6 1013 16.921 < 2.2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1#Age in years as factor with reference year being 4 (Age at t-1)
contrasts(redData$R_age_y) <- contr.treatment(levels(redData$R_age_y), base=which(levels(redData$R_age_y) == '4'))
mod1f<-lmer(as.numeric(paste(redData$diff_rltl_res))~ as.factor(redData$R_age_y) + (1|redData$recoded_id), na.action=na.exclude)## boundary (singular) fit: see ?isSingularsummary(mod1f)## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula:
## as.numeric(paste(redData$diff_rltl_res)) ~ as.factor(redData$R_age_y) +
## (1 | redData$recoded_id)
##
## REML criterion at convergence: -691.2
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -3.6127 -0.6744 0.0379 0.6816 3.5134
##
## Random effects:
## Groups Name Variance Std.Dev.
## redData$recoded_id (Intercept) 0.00000 0.0000
## Residual 0.02871 0.1694
## Number of obs: 1020, groups: redData$recoded_id, 305
##
## Fixed effects:
## Estimate Std. Error df t value Pr(>|t|)
## (Intercept) -3.707e-03 1.412e-02 1.013e+03 -0.263 0.793
## as.factor(redData$R_age_y)0 -9.686e-02 6.558e-02 1.013e+03 -1.477 0.140
## as.factor(redData$R_age_y)1 -1.114e-01 1.730e-02 1.013e+03 -6.440 1.84e-10
## as.factor(redData$R_age_y)2 -5.150e-03 1.737e-02 1.013e+03 -0.297 0.767
## as.factor(redData$R_age_y)3 1.581e-02 1.807e-02 1.013e+03 0.875 0.382
## as.factor(redData$R_age_y)5 3.080e-02 2.686e-02 1.013e+03 1.147 0.252
## as.factor(redData$R_age_y)6 -2.261e-02 4.043e-02 1.013e+03 -0.559 0.576
##
## (Intercept)
## as.factor(redData$R_age_y)0
## as.factor(redData$R_age_y)1 ***
## as.factor(redData$R_age_y)2
## as.factor(redData$R_age_y)3
## as.factor(redData$R_age_y)5
## as.factor(redData$R_age_y)6
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) a.(D$R__)0 a.(D$R__)1 a.(D$R__)2 a.(D$R__)3 a.(D$R__)5
## as.(D$R__)0 -0.215
## as.(D$R__)1 -0.816 0.176
## as.(D$R__)2 -0.813 0.175 0.664
## as.(D$R__)3 -0.782 0.168 0.638 0.635
## as.(D$R__)5 -0.526 0.113 0.429 0.427 0.411
## as.(D$R__)6 -0.349 0.075 0.285 0.284 0.273 0.184
## convergence code: 0
## boundary (singular) fit: see ?isSingularanova(mod1f)## Type III Analysis of Variance Table with Satterthwaite's method
## Sum Sq Mean Sq NumDF DenDF F value Pr(>F)
## as.factor(redData$R_age_y) 2.9149 0.48582 6 1013 16.921 < 2.2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1#Age in years as factor with reference year being 5 (Age at t-1)
contrasts(redData$R_age_y) <- contr.treatment(levels(redData$R_age_y), base=which(levels(redData$R_age_y) == '5'))
mod1f<-lmer(as.numeric(paste(redData$diff_rltl_res))~ as.factor(redData$R_age_y) + (1|redData$recoded_id), na.action=na.exclude)## boundary (singular) fit: see ?isSingularsummary(mod1f)## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula:
## as.numeric(paste(redData$diff_rltl_res)) ~ as.factor(redData$R_age_y) +
## (1 | redData$recoded_id)
##
## REML criterion at convergence: -691.2
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -3.6127 -0.6744 0.0379 0.6816 3.5134
##
## Random effects:
## Groups Name Variance Std.Dev.
## redData$recoded_id (Intercept) 0.00000 0.0000
## Residual 0.02871 0.1694
## Number of obs: 1020, groups: redData$recoded_id, 305
##
## Fixed effects:
## Estimate Std. Error df t value Pr(>|t|)
## (Intercept) 0.02709 0.02285 1013.00000 1.186 0.2360
## as.factor(redData$R_age_y)0 -0.12766 0.06800 1013.00000 -1.877 0.0608
## as.factor(redData$R_age_y)1 -0.14224 0.02494 1013.00000 -5.703 1.55e-08
## as.factor(redData$R_age_y)2 -0.03595 0.02498 1013.00000 -1.439 0.1505
## as.factor(redData$R_age_y)3 -0.01499 0.02548 1013.00000 -0.588 0.5565
## as.factor(redData$R_age_y)4 -0.03080 0.02686 1013.00000 -1.147 0.2518
## as.factor(redData$R_age_y)6 -0.05341 0.04424 1013.00000 -1.207 0.2276
##
## (Intercept)
## as.factor(redData$R_age_y)0 .
## as.factor(redData$R_age_y)1 ***
## as.factor(redData$R_age_y)2
## as.factor(redData$R_age_y)3
## as.factor(redData$R_age_y)4
## as.factor(redData$R_age_y)6
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) a.(D$R__)0 a.(D$R__)1 a.(D$R__)2 a.(D$R__)3 a.(D$R__)4
## as.(D$R__)0 -0.336
## as.(D$R__)1 -0.916 0.308
## as.(D$R__)2 -0.914 0.307 0.838
## as.(D$R__)3 -0.897 0.301 0.822 0.820
## as.(D$R__)4 -0.851 0.286 0.779 0.778 0.763
## as.(D$R__)6 -0.516 0.174 0.473 0.472 0.463 0.439
## convergence code: 0
## boundary (singular) fit: see ?isSingularanova(mod1f)## Type III Analysis of Variance Table with Satterthwaite's method
## Sum Sq Mean Sq NumDF DenDF F value Pr(>F)
## as.factor(redData$R_age_y) 2.9149 0.48582 6 1013 16.921 < 2.2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1lm<-lm(as.numeric(paste(redData$diff_rltl_res))~as.factor(redData$age_y))
summary(lm)##
## Call:
## lm(formula = as.numeric(paste(redData$diff_rltl_res)) ~ as.factor(redData$age_y))
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.61215 -0.11427 0.00643 0.11550 0.59533
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -0.10057 0.06404 -1.570 0.1167
## as.factor(redData$age_y)1 -0.01458 0.06482 -0.225 0.8221
## as.factor(redData$age_y)2 0.09171 0.06484 1.414 0.1575
## as.factor(redData$age_y)3 0.11267 0.06503 1.733 0.0835 .
## as.factor(redData$age_y)4 0.09686 0.06558 1.477 0.1400
## as.factor(redData$age_y)5 0.12766 0.06800 1.877 0.0608 .
## as.factor(redData$age_y)6 0.07425 0.07441 0.998 0.3186
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.1694 on 1013 degrees of freedom
## Multiple R-squared: 0.09109, Adjusted R-squared: 0.08571
## F-statistic: 16.92 on 6 and 1013 DF, p-value: < 2.2e-16anova(lm)## Analysis of Variance Table
##
## Response: as.numeric(paste(redData$diff_rltl_res))
## Df Sum Sq Mean Sq F value Pr(>F)
## as.factor(redData$age_y) 6 2.9149 0.48582 16.921 < 2.2e-16 ***
## Residuals 1013 29.0849 0.02871
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1nrow(subset(redData, redData$age_y==6))## [1] 20#lsmeans(lm,pairwise~as.factor(redData$age_y), adjust="Bonferroni")Test if RLTL depends on RLTL measurement at t-1
mod1<-lmer(as.numeric(paste(redData$residual_rltl_t))~ as.numeric(paste(redData$residual_rltl_tminus1))+redData$age_y + redData$feed_group + redData$genetic_group+ redData$birth_year+ as.numeric(paste(redData$sample_interval))+ redData$health_event_2week+ (1|redData$recoded_id), na.action=na.exclude)
summary(mod1)## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula:
## as.numeric(paste(redData$residual_rltl_t)) ~ as.numeric(paste(redData$residual_rltl_tminus1)) +
## redData$age_y + redData$feed_group + redData$genetic_group +
## redData$birth_year + as.numeric(paste(redData$sample_interval)) +
## redData$health_event_2week + (1 | redData$recoded_id)
##
## REML criterion at convergence: -1064.1
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -2.51548 -0.66303 -0.06712 0.55208 3.08068
##
## Random effects:
## Groups Name Variance Std.Dev.
## redData$recoded_id (Intercept) 0.005512 0.07425
## Residual 0.015346 0.12388
## Number of obs: 1020, groups: redData$recoded_id, 305
##
## Fixed effects:
## Estimate Std. Error
## (Intercept) 1.588e+01 1.126e+01
## as.numeric(paste(redData$residual_rltl_tminus1)) 1.577e-01 2.921e-02
## redData$age_y 6.145e-03 3.474e-03
## redData$feed_group2 6.969e-03 1.336e-02
## redData$feed_groupNULL 8.637e-03 1.640e-02
## redData$genetic_groupS 9.229e-03 1.239e-02
## redData$birth_year -7.928e-03 5.603e-03
## as.numeric(paste(redData$sample_interval)) 1.349e-05 4.020e-05
## redData$health_event_2week -4.652e-03 1.920e-02
## df t value Pr(>|t|)
## (Intercept) 2.079e+02 1.410 0.1601
## as.numeric(paste(redData$residual_rltl_tminus1)) 9.784e+02 5.397 8.49e-08
## redData$age_y 6.804e+02 1.769 0.0773
## redData$feed_group2 1.446e+02 0.522 0.6027
## redData$feed_groupNULL 1.787e+02 0.527 0.5990
## redData$genetic_groupS 1.531e+02 0.745 0.4573
## redData$birth_year 2.078e+02 -1.415 0.1586
## as.numeric(paste(redData$sample_interval)) 8.790e+02 0.335 0.7373
## redData$health_event_2week 8.576e+02 -0.242 0.8086
##
## (Intercept)
## as.numeric(paste(redData$residual_rltl_tminus1)) ***
## redData$age_y .
## redData$feed_group2
## redData$feed_groupNULL
## redData$genetic_groupS
## redData$birth_year
## as.numeric(paste(redData$sample_interval))
## redData$health_event_2week
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) a.((D$__ rdDt$g_ rdD$_2 rD$_NU rdD$_S rdDt$b_ a.((D$_)
## a.((D$__1)) -0.099
## redData$g_y -0.223 0.268
## rdDt$fd_gr2 -0.015 -0.037 -0.068
## rdDt$f_NULL 0.152 -0.055 0.006 0.392
## rdDt$gntc_S 0.115 -0.074 -0.012 -0.025 0.239
## rdDt$brth_y -1.000 0.099 0.223 0.015 -0.153 -0.115
## as.n((D$_)) -0.006 -0.010 -0.184 0.040 -0.027 0.053 0.004
## rdDt$hlt__2 -0.003 0.045 -0.129 0.091 0.087 0.002 0.003 0.093anova(mod1)## Type III Analysis of Variance Table with Satterthwaite's method
## Sum Sq Mean Sq NumDF DenDF
## as.numeric(paste(redData$residual_rltl_tminus1)) 0.44706 0.44706 1 978.44
## redData$age_y 0.04802 0.04802 1 680.44
## redData$feed_group 0.00606 0.00303 2 161.75
## redData$genetic_group 0.00852 0.00852 1 153.06
## redData$birth_year 0.03073 0.03073 1 207.75
## as.numeric(paste(redData$sample_interval)) 0.00173 0.00173 1 878.98
## redData$health_event_2week 0.00090 0.00090 1 857.55
## F value Pr(>F)
## as.numeric(paste(redData$residual_rltl_tminus1)) 29.1317 8.486e-08 ***
## redData$age_y 3.1294 0.07734 .
## redData$feed_group 0.1975 0.82099
## redData$genetic_group 0.5553 0.45732
## redData$birth_year 2.0022 0.15857
## as.numeric(paste(redData$sample_interval)) 0.1126 0.73733
## redData$health_event_2week 0.0587 0.80862
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1mod1<-lmer(as.numeric(paste(redData$residual_rltl_t))~ as.numeric(paste(redData$residual_rltl_tminus1))+redData$age_y + (1|redData$recoded_id), na.action=na.exclude)
summary(mod1)## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula:
## as.numeric(paste(redData$residual_rltl_t)) ~ as.numeric(paste(redData$residual_rltl_tminus1)) +
## redData$age_y + (1 | redData$recoded_id)
##
## REML criterion at convergence: -1114.7
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -2.47562 -0.65780 -0.07042 0.55074 3.06859
##
## Random effects:
## Groups Name Variance Std.Dev.
## redData$recoded_id (Intercept) 0.005188 0.07203
## Residual 0.015458 0.12433
## Number of obs: 1020, groups: redData$recoded_id, 305
##
## Fixed effects:
## Estimate Std. Error
## (Intercept) -0.047332 0.009727
## as.numeric(paste(redData$residual_rltl_tminus1)) 0.172457 0.028866
## redData$age_y 0.007691 0.003290
## df t value Pr(>|t|)
## (Intercept) 887.019161 -4.866 1.35e-06
## as.numeric(paste(redData$residual_rltl_tminus1)) 958.118312 5.974 3.25e-09
## redData$age_y 797.726818 2.337 0.0197
##
## (Intercept) ***
## as.numeric(paste(redData$residual_rltl_tminus1)) ***
## redData$age_y *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) a.((D$
## a.((D$__1)) -0.233
## redData$g_y -0.804 0.264anova(mod1)## Type III Analysis of Variance Table with Satterthwaite's method
## Sum Sq Mean Sq NumDF DenDF
## as.numeric(paste(redData$residual_rltl_tminus1)) 0.55177 0.55177 1 958.12
## redData$age_y 0.08446 0.08446 1 797.73
## F value Pr(>F)
## as.numeric(paste(redData$residual_rltl_tminus1)) 35.6935 3.253e-09 ***
## redData$age_y 5.4638 0.01966 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1mod1out<-as.data.frame(VarCorr(mod1),comp=c("Variance", "Std.Dev."))
repET<-mod1out[1,4]/sum(mod1out[1:nrow(mod1out),4])
repET## [1] 0.2512912#ggplot(redData, aes(x=age_d, y=as.numeric(paste(redData$diff_rltl_res))))+geom_point()+theme_gray(20)+ xlab("Age in days")+ ylab("RLTL change")
agePlot<- ggplot(redData, aes(x=as.factor(redData$age_y), y=as.numeric(paste(redData$diff_rltl_res))))+geom_jitter()+theme_gray(20)+ xlab("Age in years")+ ylab("RLTL change")+ geom_boxplot()
grid.arrange(SubsPlot, agePlot, ncol=2)## Warning: Use of `redData$age_y` is discouraged. Use `age_y` instead.## Warning: Use of `redData$diff_rltl_res` is discouraged. Use `diff_rltl_res`
## instead.## Warning: Use of `redData$age_y` is discouraged. Use `age_y` instead.## Warning: Use of `redData$diff_rltl_res` is discouraged. Use `diff_rltl_res`
## instead.
Investigate the effect of average telomere length over life on productive life span:
produc <-as.numeric(paste(uDeadData$amp_305_lactation))## Warning: NAs introduced by coercionmod2<-glm(as.numeric(paste(uDeadData$herd_life))~uDeadData$average_resid + produc)
summary(mod2)##
## Call:
## glm(formula = as.numeric(paste(uDeadData$herd_life)) ~ uDeadData$average_resid +
## produc)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -794.08 -361.00 -35.95 311.94 1283.75
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 1205.21924 76.16160 15.825 < 2e-16 ***
## uDeadData$average_resid -135.44248 299.38143 -0.452 0.651
## produc 0.07095 0.01075 6.597 4.18e-10 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for gaussian family taken to be 215047.8)
##
## Null deviance: 49655637 on 189 degrees of freedom
## Residual deviance: 40213930 on 187 degrees of freedom
## (51 observations deleted due to missingness)
## AIC: 2877.1
##
## Number of Fisher Scoring iterations: 2Average telomere length over life is not significantly correlated with productive lifespan.
Average telomere change over life is associated with productive lifespan.
Repeat analysis for animals with more than two samples: uDeadDataMore2
redMore2<-read.csv("/Users/luise/Documents/Work/CattleTelomereLength/LongTelomereChangePaper/git/TelomereChangeInDairyCattle/more_2_sampl.csv")
u_red_more_2<-redMore2[!duplicated(redMore2$recoded_id),]
redDeadMore2<-subset(redMore2, redMore2$herd_life != "NULL")
uDeadDataMore2<-redDeadMore2[!duplicated(redDeadMore2$recoded_id),]
product<-as.numeric(paste(uDeadDataMore2$amp_305_lactation))## Warning: NAs introduced by coercionmod2<-glm(as.numeric(paste(uDeadDataMore2$herd_life))~uDeadDataMore2$average_resid+product)
summary(mod2)##
## Call:
## glm(formula = as.numeric(paste(uDeadDataMore2$herd_life)) ~ uDeadDataMore2$average_resid +
## product)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -774.53 -377.92 -42.55 335.50 1238.72
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 1.271e+03 8.264e+01 15.382 < 2e-16 ***
## uDeadDataMore2$average_resid 3.831e+00 3.035e+02 0.013 0.99
## product 6.308e-02 1.133e-02 5.566 9.25e-08 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for gaussian family taken to be 215215.5)
##
## Null deviance: 45657680 on 183 degrees of freedom
## Residual deviance: 38954006 on 181 degrees of freedom
## (29 observations deleted due to missingness)
## AIC: 2786.6
##
## Number of Fisher Scoring iterations: 2mod3a<-glm(as.numeric(paste(uDeadDataMore2$herd_life))~ as.numeric(paste(uDeadDataMore2$rltl_change))+product)
summary(mod3a)##
## Call:
## glm(formula = as.numeric(paste(uDeadDataMore2$herd_life)) ~ as.numeric(paste(uDeadDataMore2$rltl_change)) +
## product)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -787.33 -364.01 -52.38 318.81 1224.01
##
## Coefficients:
## Estimate Std. Error t value
## (Intercept) 1.294e+03 8.446e+01 15.316
## as.numeric(paste(uDeadDataMore2$rltl_change)) 5.775e+02 6.028e+02 0.958
## product 6.238e-02 1.130e-02 5.521
## Pr(>|t|)
## (Intercept) < 2e-16 ***
## as.numeric(paste(uDeadDataMore2$rltl_change)) 0.339
## product 1.16e-07 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for gaussian family taken to be 214130)
##
## Null deviance: 45657680 on 183 degrees of freedom
## Residual deviance: 38757530 on 181 degrees of freedom
## (29 observations deleted due to missingness)
## AIC: 2785.6
##
## Number of Fisher Scoring iterations: 2AIC(mod3a)## [1] 2785.623mod3b<-glm(as.numeric(paste(uDeadDataMore2$herd_life))~ as.numeric(paste(uDeadDataMore2$abs_rltl_change))+ product)
summary(mod3b) ##
## Call:
## glm(formula = as.numeric(paste(uDeadDataMore2$herd_life)) ~ as.numeric(paste(uDeadDataMore2$abs_rltl_change)) +
## product)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -810.56 -386.38 -32.43 326.01 1251.06
##
## Coefficients:
## Estimate Std. Error t value
## (Intercept) 1344.28343 110.91823 12.120
## as.numeric(paste(uDeadDataMore2$abs_rltl_change)) -541.34365 558.66779 -0.969
## product 0.06374 0.01129 5.643
## Pr(>|t|)
## (Intercept) < 2e-16 ***
## as.numeric(paste(uDeadDataMore2$abs_rltl_change)) 0.334
## product 6.34e-08 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for gaussian family taken to be 214105)
##
## Null deviance: 45657680 on 183 degrees of freedom
## Residual deviance: 38753008 on 181 degrees of freedom
## (29 observations deleted due to missingness)
## AIC: 2785.6
##
## Number of Fisher Scoring iterations: 2AIC(mod3b)## [1] 2785.601AIC(mod3b)-AIC(mod3a)## [1] -0.02147028#All together
mod3c<-glm(as.numeric(paste(uDeadDataMore2$herd_life))~ as.numeric(paste(uDeadDataMore2$rltl_change))+ as.numeric(paste(uDeadDataMore2$abs_rltl_change))+ uDeadDataMore2$average_resid +product)
summary(mod3c)##
## Call:
## glm(formula = as.numeric(paste(uDeadDataMore2$herd_life)) ~ as.numeric(paste(uDeadDataMore2$rltl_change)) +
## as.numeric(paste(uDeadDataMore2$abs_rltl_change)) + uDeadDataMore2$average_resid +
## product)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -813.50 -364.81 -61.04 328.87 1236.41
##
## Coefficients:
## Estimate Std. Error t value
## (Intercept) 1345.44375 114.08909 11.793
## as.numeric(paste(uDeadDataMore2$rltl_change)) 443.74881 644.06446 0.689
## as.numeric(paste(uDeadDataMore2$abs_rltl_change)) -421.00056 595.89616 -0.706
## uDeadDataMore2$average_resid 1.53645 311.59944 0.005
## product 0.06304 0.01140 5.532
## Pr(>|t|)
## (Intercept) < 2e-16 ***
## as.numeric(paste(uDeadDataMore2$rltl_change)) 0.492
## as.numeric(paste(uDeadDataMore2$abs_rltl_change)) 0.481
## uDeadDataMore2$average_resid 0.996
## product 1.11e-07 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for gaussian family taken to be 215903.2)
##
## Null deviance: 45657680 on 183 degrees of freedom
## Residual deviance: 38646678 on 179 degrees of freedom
## (29 observations deleted due to missingness)
## AIC: 2789.1
##
## Number of Fisher Scoring iterations: 2mod3d<-glm(as.numeric(paste(uDeadDataMore2$herd_life))~ as.numeric(paste(uDeadDataMore2$rltl_change))+ as.numeric(paste(uDeadDataMore2$abs_rltl_change))+ product)
summary(mod3d)##
## Call:
## glm(formula = as.numeric(paste(uDeadDataMore2$herd_life)) ~ as.numeric(paste(uDeadDataMore2$rltl_change)) +
## as.numeric(paste(uDeadDataMore2$abs_rltl_change)) + product)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -813.44 -364.72 -60.94 328.88 1236.64
##
## Coefficients:
## Estimate Std. Error t value
## (Intercept) 1345.32218 111.08303 12.111
## as.numeric(paste(uDeadDataMore2$rltl_change)) 444.33049 631.40750 0.704
## as.numeric(paste(uDeadDataMore2$abs_rltl_change)) -420.49044 585.21439 -0.719
## product 0.06305 0.01135 5.554
## Pr(>|t|)
## (Intercept) < 2e-16 ***
## as.numeric(paste(uDeadDataMore2$rltl_change)) 0.483
## as.numeric(paste(uDeadDataMore2$abs_rltl_change)) 0.473
## product 9.92e-08 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for gaussian family taken to be 214703.8)
##
## Null deviance: 45657680 on 183 degrees of freedom
## Residual deviance: 38646683 on 180 degrees of freedom
## (29 observations deleted due to missingness)
## AIC: 2787.1
##
## Number of Fisher Scoring iterations: 2Cox proportional hazard model of productive lifespan with average lifetime RLTL as predictor
coxFit<-coxph(Surv(coxData$TimeMeasure, coxData$Event == 1)~coxData$average_resid)
summary(coxFit)## Call:
## coxph(formula = Surv(coxData$TimeMeasure, coxData$Event == 1) ~
## coxData$average_resid)
##
## n= 305, number of events= 241
##
## coef exp(coef) se(coef) z Pr(>|z|)
## coxData$average_resid 0.3407 1.4059 0.5905 0.577 0.564
##
## exp(coef) exp(-coef) lower .95 upper .95
## coxData$average_resid 1.406 0.7113 0.4419 4.473
##
## Concordance= 0.509 (se = 0.021 )
## Likelihood ratio test= 0.33 on 1 df, p=0.6
## Wald test = 0.33 on 1 df, p=0.6
## Score (logrank) test = 0.33 on 1 df, p=0.6summary(coxFit)$logtest["pvalue"]## pvalue
## 0.5641153Change to discrete scale to enable plotting
### discrete scale
coxFit<-coxph(Surv(coxData$TimeMeasure, coxData$Event == 1)~coxData$average_resid_3G)
summary(coxFit)## Call:
## coxph(formula = Surv(coxData$TimeMeasure, coxData$Event == 1) ~
## coxData$average_resid_3G)
##
## n= 305, number of events= 241
##
## coef exp(coef) se(coef) z Pr(>|z|)
## coxData$average_resid_3G 0.01410 1.01420 0.07886 0.179 0.858
##
## exp(coef) exp(-coef) lower .95 upper .95
## coxData$average_resid_3G 1.014 0.986 0.869 1.184
##
## Concordance= 0.501 (se = 0.02 )
## Likelihood ratio test= 0.03 on 1 df, p=0.9
## Wald test = 0.03 on 1 df, p=0.9
## Score (logrank) test = 0.03 on 1 df, p=0.9#ggsave(file = "/Users/luise/Documents/Work/CattleTelomereLength/AgingCell/Images/MeanRLTL.pdf", print(g2))
#write.csv(coxData, "/Users/luise/Documents/Work/CattleTelomereLength/AgingCell/rawdata/CaxData.csv")Simple correlations of herd life with lifetime RLTL measures
They show the same as Cox proportional hazard models.
dData<-subset(data, data$herd_life != "NULL")
cor.test(as.numeric(paste(dData$herd_life)), dData$abs_rltl_change)##
## Pearson's product-moment correlation
##
## data: as.numeric(paste(dData$herd_life)) and dData$abs_rltl_change
## t = -5.5853, df = 1020, p-value = 2.991e-08
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
## -0.2311491 -0.1121316
## sample estimates:
## cor
## -0.172269cor.test(as.numeric(paste(dData$herd_life)), dData$rltl_change)##
## Pearson's product-moment correlation
##
## data: as.numeric(paste(dData$herd_life)) and dData$rltl_change
## t = 8.4519, df = 1020, p-value < 2.2e-16
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
## 0.1976120 0.3122569
## sample estimates:
## cor
## 0.2558338cor.test(as.numeric(paste(dData$herd_life)), dData$average_resid)##
## Pearson's product-moment correlation
##
## data: as.numeric(paste(dData$herd_life)) and dData$average_resid
## t = -0.19531, df = 1020, p-value = 0.8452
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
## -0.06741195 0.05522731
## sample estimates:
## cor
## -0.006115315### discrete scale
coxFit<-coxph(Surv(coxData$TimeMeasure, coxData$Event == 1)~coxData$abs_rltl_change_3G)
summary(coxFit)## Call:
## coxph(formula = Surv(coxData$TimeMeasure, coxData$Event == 1) ~
## coxData$abs_rltl_change_3G)
##
## n= 305, number of events= 241
##
## coef exp(coef) se(coef) z Pr(>|z|)
## coxData$abs_rltl_change_3G 0.17936 1.19645 0.08189 2.19 0.0285 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## exp(coef) exp(-coef) lower .95 upper .95
## coxData$abs_rltl_change_3G 1.196 0.8358 1.019 1.405
##
## Concordance= 0.542 (se = 0.02 )
## Likelihood ratio test= 4.8 on 1 df, p=0.03
## Wald test = 4.8 on 1 df, p=0.03
## Score (logrank) test = 4.82 on 1 df, p=0.03coxFit<-survfit(Surv(coxData$TimeMeasure, coxData$Event == 1)~coxData$abs_rltl_change_3G)
summary(coxFit)## Call: survfit(formula = Surv(coxData$TimeMeasure, coxData$Event ==
## 1) ~ coxData$abs_rltl_change_3G)
##
## coxData$abs_rltl_change_3G=1
## time n.risk n.event survival std.err lower 95% CI upper 95% CI
## 495 102 1 0.990 0.00976 0.9713 1.000
## 501 101 1 0.980 0.01373 0.9539 1.000
## 636 100 1 0.971 0.01673 0.9383 1.000
## 641 99 1 0.961 0.01922 0.9238 0.999
## 648 98 1 0.951 0.02138 0.9100 0.994
## 705 97 1 0.941 0.02330 0.8966 0.988
## 724 96 1 0.931 0.02503 0.8836 0.982
## 727 95 1 0.922 0.02662 0.8708 0.975
## 735 94 1 0.912 0.02808 0.8583 0.969
## 794 93 1 0.902 0.02944 0.8461 0.962
## 805 92 1 0.892 0.03071 0.8339 0.954
## 830 91 1 0.882 0.03190 0.8220 0.947
## 831 90 1 0.873 0.03302 0.8102 0.940
## 833 89 1 0.863 0.03407 0.7985 0.932
## 836 88 1 0.853 0.03507 0.7869 0.925
## 895 87 1 0.843 0.03601 0.7754 0.917
## 897 86 1 0.833 0.03690 0.7641 0.909
## 904 85 1 0.824 0.03775 0.7528 0.901
## 909 84 1 0.814 0.03855 0.7416 0.893
## 956 83 1 0.804 0.03931 0.7304 0.885
## 995 82 1 0.794 0.04004 0.7194 0.877
## 1033 81 1 0.784 0.04072 0.7084 0.868
## 1082 80 1 0.775 0.04138 0.6975 0.860
## 1121 79 1 0.765 0.04200 0.6867 0.852
## 1128 78 1 0.755 0.04259 0.6759 0.843
## 1137 77 1 0.745 0.04315 0.6651 0.835
## 1158 76 2 0.725 0.04419 0.6439 0.817
## 1186 74 1 0.716 0.04466 0.6333 0.809
## 1209 73 1 0.706 0.04512 0.6228 0.800
## 1341 72 1 0.696 0.04554 0.6123 0.791
## 1356 71 1 0.686 0.04594 0.6019 0.782
## 1381 70 1 0.676 0.04632 0.5915 0.774
## 1483 69 1 0.667 0.04668 0.5812 0.765
## 1515 68 1 0.657 0.04701 0.5709 0.756
## 1551 67 1 0.647 0.04732 0.5607 0.747
## 1567 66 1 0.637 0.04761 0.5505 0.738
## 1573 65 1 0.627 0.04787 0.5403 0.729
## 1588 64 1 0.618 0.04812 0.5302 0.720
## 1596 63 1 0.608 0.04834 0.5201 0.710
## 1604 62 1 0.598 0.04855 0.5101 0.701
## 1612 61 1 0.588 0.04873 0.5001 0.692
## 1620 60 1 0.578 0.04889 0.4901 0.683
## 1624 59 1 0.569 0.04904 0.4802 0.673
## 1627 58 1 0.559 0.04916 0.4703 0.664
## 1633 57 1 0.549 0.04927 0.4605 0.655
## 1657 56 1 0.539 0.04935 0.4507 0.645
## 1766 55 1 0.529 0.04942 0.4409 0.636
## 1815 54 1 0.520 0.04947 0.4312 0.626
## 1825 53 1 0.510 0.04950 0.4215 0.617
## 1830 52 1 0.500 0.04951 0.4118 0.607
## 1834 51 1 0.490 0.04950 0.4022 0.597
## 1902 50 1 0.480 0.04947 0.3926 0.588
## 1903 49 1 0.471 0.04942 0.3830 0.578
## 1942 48 1 0.461 0.04935 0.3735 0.568
## 1960 47 1 0.451 0.04927 0.3641 0.559
## 2015 43 1 0.440 0.04923 0.3538 0.548
## 2017 42 1 0.430 0.04916 0.3437 0.538
## 2060 40 1 0.419 0.04909 0.3333 0.527
## 2072 39 1 0.409 0.04900 0.3229 0.517
## 2077 37 1 0.397 0.04890 0.3123 0.506
## 2140 35 1 0.386 0.04880 0.3014 0.495
## 2146 33 1 0.374 0.04871 0.2901 0.483
## 2177 31 1 0.362 0.04861 0.2786 0.471
## 2192 29 1 0.350 0.04851 0.2666 0.459
## 2214 28 1 0.337 0.04836 0.2547 0.447
## 2264 27 1 0.325 0.04816 0.2429 0.434
## 2275 26 1 0.312 0.04790 0.2313 0.422
## 2375 19 1 0.296 0.04812 0.2152 0.407
## 2437 15 1 0.276 0.04879 0.1954 0.390
## 2452 13 1 0.255 0.04944 0.1743 0.373
## 2457 12 1 0.234 0.04968 0.1541 0.354
## 2476 11 1 0.212 0.04950 0.1346 0.335
## 2567 7 1 0.182 0.05089 0.1053 0.315
## 2627 5 1 0.146 0.05214 0.0722 0.294
## 2708 4 1 0.109 0.05024 0.0444 0.269
##
## coxData$abs_rltl_change_3G=2
## time n.risk n.event survival std.err lower 95% CI upper 95% CI
## 383 102 1 0.9902 0.00976 0.9713 1.000
## 499 101 1 0.9804 0.01373 0.9539 1.000
## 628 100 1 0.9706 0.01673 0.9383 1.000
## 658 99 1 0.9608 0.01922 0.9238 0.999
## 721 98 1 0.9510 0.02138 0.9100 0.994
## 731 97 1 0.9412 0.02330 0.8966 0.988
## 749 96 1 0.9314 0.02503 0.8836 0.982
## 822 95 1 0.9216 0.02662 0.8708 0.975
## 840 94 1 0.9118 0.02808 0.8583 0.969
## 908 93 1 0.9020 0.02944 0.8461 0.962
## 923 92 1 0.8922 0.03071 0.8339 0.954
## 928 91 1 0.8824 0.03190 0.8220 0.947
## 930 90 1 0.8725 0.03302 0.8102 0.940
## 948 89 1 0.8627 0.03407 0.7985 0.932
## 961 88 1 0.8529 0.03507 0.7869 0.925
## 1019 87 1 0.8431 0.03601 0.7754 0.917
## 1057 86 1 0.8333 0.03690 0.7641 0.909
## 1094 85 1 0.8235 0.03775 0.7528 0.901
## 1110 84 1 0.8137 0.03855 0.7416 0.893
## 1127 83 1 0.8039 0.03931 0.7304 0.885
## 1138 82 1 0.7941 0.04004 0.7194 0.877
## 1142 81 1 0.7843 0.04072 0.7084 0.868
## 1160 80 1 0.7745 0.04138 0.6975 0.860
## 1182 79 1 0.7647 0.04200 0.6867 0.852
## 1193 78 1 0.7549 0.04259 0.6759 0.843
## 1263 77 1 0.7451 0.04315 0.6651 0.835
## 1280 76 1 0.7353 0.04368 0.6545 0.826
## 1283 75 1 0.7255 0.04419 0.6439 0.817
## 1371 74 1 0.7157 0.04466 0.6333 0.809
## 1417 73 1 0.7059 0.04512 0.6228 0.800
## 1421 72 1 0.6961 0.04554 0.6123 0.791
## 1429 71 1 0.6863 0.04594 0.6019 0.782
## 1461 70 1 0.6765 0.04632 0.5915 0.774
## 1468 69 1 0.6667 0.04668 0.5812 0.765
## 1479 68 1 0.6569 0.04701 0.5709 0.756
## 1516 67 1 0.6471 0.04732 0.5607 0.747
## 1532 66 1 0.6373 0.04761 0.5505 0.738
## 1549 65 1 0.6275 0.04787 0.5403 0.729
## 1575 64 1 0.6176 0.04812 0.5302 0.720
## 1593 63 1 0.6078 0.04834 0.5201 0.710
## 1612 62 1 0.5980 0.04855 0.5101 0.701
## 1633 61 1 0.5882 0.04873 0.5001 0.692
## 1638 60 1 0.5784 0.04889 0.4901 0.683
## 1661 59 1 0.5686 0.04904 0.4802 0.673
## 1693 58 1 0.5588 0.04916 0.4703 0.664
## 1739 57 1 0.5490 0.04927 0.4605 0.655
## 1750 56 1 0.5392 0.04935 0.4507 0.645
## 1758 55 1 0.5294 0.04942 0.4409 0.636
## 1776 54 1 0.5196 0.04947 0.4312 0.626
## 1785 53 1 0.5098 0.04950 0.4215 0.617
## 1786 52 1 0.5000 0.04951 0.4118 0.607
## 1790 51 1 0.4902 0.04950 0.4022 0.597
## 1793 50 1 0.4804 0.04947 0.3926 0.588
## 1822 49 1 0.4706 0.04942 0.3830 0.578
## 1855 48 1 0.4608 0.04935 0.3735 0.568
## 1904 47 1 0.4510 0.04927 0.3641 0.559
## 1914 46 1 0.4412 0.04916 0.3546 0.549
## 1927 45 1 0.4314 0.04904 0.3452 0.539
## 1975 43 1 0.4213 0.04891 0.3356 0.529
## 1977 42 1 0.4113 0.04877 0.3260 0.519
## 1978 41 1 0.4013 0.04860 0.3165 0.509
## 1992 40 1 0.3912 0.04841 0.3070 0.499
## 2018 39 1 0.3812 0.04819 0.2975 0.488
## 2041 38 1 0.3712 0.04796 0.2881 0.478
## 2045 37 1 0.3611 0.04770 0.2788 0.468
## 2063 36 1 0.3511 0.04742 0.2695 0.458
## 2079 35 1 0.3411 0.04711 0.2602 0.447
## 2090 34 1 0.3311 0.04678 0.2510 0.437
## 2155 28 1 0.3192 0.04658 0.2398 0.425
## 2171 27 1 0.3074 0.04633 0.2288 0.413
## 2185 26 1 0.2956 0.04604 0.2178 0.401
## 2191 25 1 0.2838 0.04569 0.2070 0.389
## 2275 22 1 0.2709 0.04539 0.1950 0.376
## 2325 20 1 0.2573 0.04510 0.1825 0.363
## 2510 17 1 0.2422 0.04492 0.1684 0.348
## 2541 16 1 0.2270 0.04459 0.1545 0.334
## 2583 13 1 0.2096 0.04445 0.1383 0.318
## 2623 11 1 0.1905 0.04430 0.1208 0.301
## 2766 8 1 0.1667 0.04471 0.0986 0.282
## 2767 7 1 0.1429 0.04421 0.0779 0.262
## 2822 5 1 0.1143 0.04364 0.0541 0.242
## 2823 4 1 0.0857 0.04103 0.0336 0.219
##
## coxData$abs_rltl_change_3G=3
## time n.risk n.event survival std.err lower 95% CI upper 95% CI
## 30 101 1 0.9901 0.00985 0.9710 1.000
## 273 100 1 0.9802 0.01386 0.9534 1.000
## 302 99 1 0.9703 0.01689 0.9377 1.000
## 330 98 1 0.9604 0.01941 0.9231 0.999
## 349 97 1 0.9505 0.02158 0.9091 0.994
## 476 96 1 0.9406 0.02352 0.8956 0.988
## 535 95 1 0.9307 0.02527 0.8825 0.982
## 544 94 1 0.9208 0.02687 0.8696 0.975
## 658 93 1 0.9109 0.02835 0.8570 0.968
## 699 92 1 0.9010 0.02972 0.8446 0.961
## 742 91 1 0.8911 0.03100 0.8324 0.954
## 750 90 1 0.8812 0.03220 0.8203 0.947
## 759 89 1 0.8713 0.03332 0.8084 0.939
## 760 88 1 0.8614 0.03438 0.7966 0.931
## 775 87 1 0.8515 0.03538 0.7849 0.924
## 798 86 1 0.8416 0.03633 0.7733 0.916
## 803 85 1 0.8317 0.03723 0.7618 0.908
## 845 84 1 0.8218 0.03808 0.7504 0.900
## 867 83 1 0.8119 0.03889 0.7391 0.892
## 888 82 1 0.8020 0.03965 0.7279 0.884
## 906 81 1 0.7921 0.04038 0.7168 0.875
## 919 80 1 0.7822 0.04107 0.7057 0.867
## 925 79 1 0.7723 0.04173 0.6947 0.859
## 932 78 1 0.7624 0.04235 0.6837 0.850
## 944 77 1 0.7525 0.04294 0.6728 0.842
## 958 76 1 0.7426 0.04350 0.6620 0.833
## 961 75 1 0.7327 0.04404 0.6513 0.824
## 982 74 1 0.7228 0.04454 0.6405 0.816
## 1006 73 1 0.7129 0.04502 0.6299 0.807
## 1045 72 1 0.7030 0.04547 0.6193 0.798
## 1069 71 1 0.6931 0.04589 0.6087 0.789
## 1142 70 1 0.6832 0.04629 0.5982 0.780
## 1147 69 1 0.6733 0.04667 0.5877 0.771
## 1156 68 1 0.6634 0.04702 0.5773 0.762
## 1165 67 1 0.6535 0.04735 0.5669 0.753
## 1177 66 1 0.6436 0.04766 0.5566 0.744
## 1195 65 1 0.6337 0.04794 0.5463 0.735
## 1205 64 1 0.6238 0.04820 0.5361 0.726
## 1222 63 1 0.6139 0.04844 0.5259 0.717
## 1301 62 1 0.6040 0.04866 0.5157 0.707
## 1345 61 1 0.5941 0.04886 0.5056 0.698
## 1364 60 1 0.5842 0.04904 0.4955 0.689
## 1373 59 1 0.5743 0.04920 0.4855 0.679
## 1378 58 1 0.5644 0.04934 0.4755 0.670
## 1384 57 1 0.5545 0.04946 0.4655 0.660
## 1436 56 1 0.5446 0.04955 0.4556 0.651
## 1530 55 1 0.5347 0.04963 0.4457 0.641
## 1557 54 1 0.5248 0.04969 0.4359 0.632
## 1561 53 1 0.5149 0.04973 0.4261 0.622
## 1565 52 1 0.5050 0.04975 0.4163 0.613
## 1597 51 1 0.4950 0.04975 0.4065 0.603
## 1608 50 1 0.4851 0.04973 0.3968 0.593
## 1609 49 1 0.4752 0.04969 0.3872 0.583
## 1610 48 1 0.4653 0.04963 0.3776 0.574
## 1615 47 1 0.4554 0.04955 0.3680 0.564
## 1671 45 1 0.4453 0.04948 0.3582 0.554
## 1701 44 1 0.4352 0.04938 0.3484 0.544
## 1715 43 1 0.4251 0.04925 0.3387 0.533
## 1756 42 1 0.4150 0.04911 0.3291 0.523
## 1762 41 1 0.4048 0.04894 0.3194 0.513
## 1766 40 1 0.3947 0.04876 0.3098 0.503
## 1770 39 1 0.3846 0.04854 0.3003 0.493
## 1786 38 1 0.3745 0.04831 0.2908 0.482
## 1848 37 1 0.3644 0.04805 0.2814 0.472
## 1904 36 1 0.3542 0.04777 0.2720 0.461
## 1937 35 1 0.3441 0.04747 0.2626 0.451
## 1940 34 1 0.3340 0.04714 0.2533 0.440
## 1994 33 1 0.3239 0.04678 0.2440 0.430
## 1999 32 1 0.3138 0.04640 0.2348 0.419
## 2001 31 1 0.3036 0.04600 0.2256 0.409
## 2035 27 1 0.2924 0.04565 0.2153 0.397
## 2060 26 1 0.2811 0.04526 0.2051 0.385
## 2067 25 2 0.2586 0.04434 0.1848 0.362
## 2075 23 1 0.2474 0.04382 0.1748 0.350
## 2086 22 1 0.2362 0.04324 0.1649 0.338
## 2113 18 1 0.2230 0.04279 0.1531 0.325
## 2166 16 1 0.2091 0.04232 0.1406 0.311
## 2235 14 1 0.1942 0.04185 0.1273 0.296
## 2285 11 1 0.1765 0.04160 0.1112 0.280
## 2333 8 1 0.1544 0.04185 0.0908 0.263
## 2486 4 1 0.1158 0.04586 0.0533 0.252
## 2579 3 1 0.0772 0.04392 0.0253 0.235#### plot results for discrete scale
g3<- ggsurvplot(coxFit, data=coxData, pval = 0.029, ggtheme = theme_minimal(), risk.table.y.text.col = T, risk.table.y.text = FALSE, legend.labs=c("Tertile 1 (least absolute RLTL change)", "Tertile 2", "Tertile 3 (most absolute RLTL change)")) +theme_survminer(40)+ guides(colour = guide_legend(nrow = 3))
ggsurvplot(coxFit, data=coxData, pval = 0.029, ggtheme = theme_minimal(), risk.table.y.text.col = T, risk.table.y.text = FALSE, legend.labs=c("Tertile 1 (least abs RLTL change)", "Tertile 2", "Tertile 3 (most abs RLTL change)")) +theme_survminer(40)
#ggsave(file = "/Users/luise/Documents/Work/CattleTelomereLength/AgingCell/Images/AbsoluteRLTLchange.pdf", print(g3))#continous scale
coxFit<-coxph(Surv(coxData$TimeMeasure, coxData$Event == 1)~coxData$rltl_change)
summary(coxFit)## Call:
## coxph(formula = Surv(coxData$TimeMeasure, coxData$Event == 1) ~
## coxData$rltl_change)
##
## n= 305, number of events= 241
##
## coef exp(coef) se(coef) z Pr(>|z|)
## coxData$rltl_change -5.209383 0.005465 0.844821 -6.166 6.99e-10 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## exp(coef) exp(-coef) lower .95 upper .95
## coxData$rltl_change 0.005465 183 0.001043 0.02862
##
## Concordance= 0.575 (se = 0.023 )
## Likelihood ratio test= 28.66 on 1 df, p=9e-08
## Wald test = 38.02 on 1 df, p=7e-10
## Score (logrank) test = 33.59 on 1 df, p=7e-09#3 groups:
coxFit<-coxph(Surv(coxData$TimeMeasure, coxData$Event == 1)~coxData$rltl_change_3G)
summary(coxFit)## Call:
## coxph(formula = Surv(coxData$TimeMeasure, coxData$Event == 1) ~
## coxData$rltl_change_3G)
##
## n= 305, number of events= 241
##
## coef exp(coef) se(coef) z Pr(>|z|)
## coxData$rltl_change_3G -0.25745 0.77302 0.08656 -2.974 0.00294 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## exp(coef) exp(-coef) lower .95 upper .95
## coxData$rltl_change_3G 0.773 1.294 0.6524 0.9159
##
## Concordance= 0.56 (se = 0.019 )
## Likelihood ratio test= 8.88 on 1 df, p=0.003
## Wald test = 8.85 on 1 df, p=0.003
## Score (logrank) test = 8.92 on 1 df, p=0.003coxFit<-survfit(Surv(coxData$TimeMeasure, coxData$Event == 1)~coxData$rltl_change_3G)
summary(coxFit)## Call: survfit(formula = Surv(coxData$TimeMeasure, coxData$Event ==
## 1) ~ coxData$rltl_change_3G)
##
## coxData$rltl_change_3G=1
## time n.risk n.event survival std.err lower 95% CI upper 95% CI
## 30 102 1 0.9902 0.00976 0.9713 1.000
## 273 101 1 0.9804 0.01373 0.9539 1.000
## 330 100 1 0.9706 0.01673 0.9383 1.000
## 349 99 1 0.9608 0.01922 0.9238 0.999
## 383 98 1 0.9510 0.02138 0.9100 0.994
## 476 97 1 0.9412 0.02330 0.8966 0.988
## 499 96 1 0.9314 0.02503 0.8836 0.982
## 501 95 1 0.9216 0.02662 0.8708 0.975
## 535 94 1 0.9118 0.02808 0.8583 0.969
## 544 93 1 0.9020 0.02944 0.8461 0.962
## 628 92 1 0.8922 0.03071 0.8339 0.954
## 636 91 1 0.8824 0.03190 0.8220 0.947
## 658 90 1 0.8725 0.03302 0.8102 0.940
## 699 89 1 0.8627 0.03407 0.7985 0.932
## 705 88 1 0.8529 0.03507 0.7869 0.925
## 724 87 1 0.8431 0.03601 0.7754 0.917
## 742 86 1 0.8333 0.03690 0.7641 0.909
## 750 85 1 0.8235 0.03775 0.7528 0.901
## 759 84 1 0.8137 0.03855 0.7416 0.893
## 760 83 1 0.8039 0.03931 0.7304 0.885
## 775 82 1 0.7941 0.04004 0.7194 0.877
## 845 81 1 0.7843 0.04072 0.7084 0.868
## 888 80 1 0.7745 0.04138 0.6975 0.860
## 906 79 1 0.7647 0.04200 0.6867 0.852
## 908 78 1 0.7549 0.04259 0.6759 0.843
## 923 77 1 0.7451 0.04315 0.6651 0.835
## 925 76 1 0.7353 0.04368 0.6545 0.826
## 944 75 1 0.7255 0.04419 0.6439 0.817
## 948 74 1 0.7157 0.04466 0.6333 0.809
## 961 73 2 0.6961 0.04554 0.6123 0.791
## 995 71 1 0.6863 0.04594 0.6019 0.782
## 1006 70 1 0.6765 0.04632 0.5915 0.774
## 1019 69 1 0.6667 0.04668 0.5812 0.765
## 1069 68 1 0.6569 0.04701 0.5709 0.756
## 1082 67 1 0.6471 0.04732 0.5607 0.747
## 1094 66 1 0.6373 0.04761 0.5505 0.738
## 1110 65 1 0.6275 0.04787 0.5403 0.729
## 1121 64 1 0.6176 0.04812 0.5302 0.720
## 1156 63 1 0.6078 0.04834 0.5201 0.710
## 1160 62 1 0.5980 0.04855 0.5101 0.701
## 1165 61 1 0.5882 0.04873 0.5001 0.692
## 1177 60 1 0.5784 0.04889 0.4901 0.683
## 1182 59 1 0.5686 0.04904 0.4802 0.673
## 1195 58 1 0.5588 0.04916 0.4703 0.664
## 1205 57 1 0.5490 0.04927 0.4605 0.655
## 1280 56 1 0.5392 0.04935 0.4507 0.645
## 1345 55 1 0.5294 0.04942 0.4409 0.636
## 1364 54 1 0.5196 0.04947 0.4312 0.626
## 1371 53 1 0.5098 0.04950 0.4215 0.617
## 1384 52 1 0.5000 0.04951 0.4118 0.607
## 1417 51 1 0.4902 0.04950 0.4022 0.597
## 1421 50 1 0.4804 0.04947 0.3926 0.588
## 1436 49 1 0.4706 0.04942 0.3830 0.578
## 1530 48 1 0.4608 0.04935 0.3735 0.568
## 1549 47 1 0.4510 0.04927 0.3641 0.559
## 1565 46 1 0.4412 0.04916 0.3546 0.549
## 1567 45 1 0.4314 0.04904 0.3452 0.539
## 1573 44 1 0.4216 0.04889 0.3358 0.529
## 1588 43 1 0.4118 0.04873 0.3265 0.519
## 1597 42 1 0.4020 0.04855 0.3172 0.509
## 1608 41 1 0.3922 0.04834 0.3080 0.499
## 1609 40 1 0.3824 0.04812 0.2988 0.489
## 1615 39 1 0.3725 0.04787 0.2896 0.479
## 1657 38 1 0.3627 0.04761 0.2805 0.469
## 1693 37 1 0.3529 0.04732 0.2714 0.459
## 1739 36 1 0.3431 0.04701 0.2623 0.449
## 1750 35 1 0.3333 0.04668 0.2533 0.439
## 1756 34 1 0.3235 0.04632 0.2444 0.428
## 1758 33 1 0.3137 0.04594 0.2354 0.418
## 1762 32 1 0.3039 0.04554 0.2266 0.408
## 1766 31 1 0.2941 0.04512 0.2177 0.397
## 1822 30 1 0.2843 0.04466 0.2090 0.387
## 1848 29 1 0.2745 0.04419 0.2002 0.376
## 1903 28 1 0.2647 0.04368 0.1916 0.366
## 1904 27 1 0.2549 0.04315 0.1829 0.355
## 1927 26 1 0.2451 0.04259 0.1744 0.345
## 1940 25 1 0.2353 0.04200 0.1658 0.334
## 1999 23 1 0.2251 0.04140 0.1569 0.323
## 2001 22 1 0.2148 0.04076 0.1481 0.312
## 2041 20 1 0.2041 0.04012 0.1388 0.300
## 2067 19 2 0.1826 0.03866 0.1206 0.277
## 2072 17 1 0.1719 0.03785 0.1116 0.265
## 2086 15 1 0.1604 0.03702 0.1020 0.252
## 2113 12 1 0.1470 0.03627 0.0907 0.238
## 2235 8 1 0.1287 0.03609 0.0742 0.223
## 2623 2 1 0.0643 0.04894 0.0145 0.286
## 2627 1 1 0.0000 NaN NA NA
##
## coxData$rltl_change_3G=2
## time n.risk n.event survival std.err lower 95% CI upper 95% CI
## 641 102 1 0.990 0.00976 0.9713 1.000
## 727 101 1 0.980 0.01373 0.9539 1.000
## 731 100 1 0.971 0.01673 0.9383 1.000
## 749 99 1 0.961 0.01922 0.9238 0.999
## 794 98 1 0.951 0.02138 0.9100 0.994
## 805 97 1 0.941 0.02330 0.8966 0.988
## 904 96 1 0.931 0.02503 0.8836 0.982
## 909 95 1 0.922 0.02662 0.8708 0.975
## 928 94 1 0.912 0.02808 0.8583 0.969
## 956 93 1 0.902 0.02944 0.8461 0.962
## 958 92 1 0.892 0.03071 0.8339 0.954
## 1127 91 1 0.882 0.03190 0.8220 0.947
## 1137 90 1 0.873 0.03302 0.8102 0.940
## 1138 89 1 0.863 0.03407 0.7985 0.932
## 1142 88 2 0.843 0.03601 0.7754 0.917
## 1158 86 1 0.833 0.03690 0.7641 0.909
## 1222 85 1 0.824 0.03775 0.7528 0.901
## 1283 84 1 0.814 0.03855 0.7416 0.893
## 1341 83 1 0.804 0.03931 0.7304 0.885
## 1356 82 1 0.794 0.04004 0.7194 0.877
## 1378 81 1 0.784 0.04072 0.7084 0.868
## 1461 80 1 0.775 0.04138 0.6975 0.860
## 1468 79 1 0.765 0.04200 0.6867 0.852
## 1479 78 1 0.755 0.04259 0.6759 0.843
## 1532 77 1 0.745 0.04315 0.6651 0.835
## 1557 76 1 0.735 0.04368 0.6545 0.826
## 1596 75 1 0.725 0.04419 0.6439 0.817
## 1610 74 1 0.716 0.04466 0.6333 0.809
## 1612 73 1 0.706 0.04512 0.6228 0.800
## 1624 72 1 0.696 0.04554 0.6123 0.791
## 1627 71 1 0.686 0.04594 0.6019 0.782
## 1633 70 1 0.676 0.04632 0.5915 0.774
## 1661 69 1 0.667 0.04668 0.5812 0.765
## 1715 68 1 0.657 0.04701 0.5709 0.756
## 1766 67 1 0.647 0.04732 0.5607 0.747
## 1776 66 1 0.637 0.04761 0.5505 0.738
## 1786 65 1 0.627 0.04787 0.5403 0.729
## 1790 64 1 0.618 0.04812 0.5302 0.720
## 1793 63 1 0.608 0.04834 0.5201 0.710
## 1815 62 1 0.598 0.04855 0.5101 0.701
## 1825 61 1 0.588 0.04873 0.5001 0.692
## 1855 60 1 0.578 0.04889 0.4901 0.683
## 1937 59 1 0.569 0.04904 0.4802 0.673
## 1942 58 1 0.559 0.04916 0.4703 0.664
## 1977 57 1 0.549 0.04927 0.4605 0.655
## 1978 56 1 0.539 0.04935 0.4507 0.645
## 1992 55 1 0.529 0.04942 0.4409 0.636
## 2015 53 1 0.519 0.04949 0.4309 0.626
## 2018 52 1 0.509 0.04953 0.4210 0.616
## 2045 51 1 0.499 0.04956 0.4112 0.607
## 2060 50 1 0.489 0.04957 0.4013 0.597
## 2075 49 1 0.479 0.04955 0.3916 0.587
## 2079 48 1 0.469 0.04951 0.3818 0.577
## 2146 44 1 0.459 0.04953 0.3713 0.567
## 2166 43 1 0.448 0.04951 0.3609 0.556
## 2171 42 1 0.437 0.04947 0.3505 0.546
## 2177 41 1 0.427 0.04940 0.3402 0.535
## 2185 40 1 0.416 0.04930 0.3299 0.525
## 2191 38 1 0.405 0.04921 0.3194 0.514
## 2264 36 1 0.394 0.04911 0.3085 0.503
## 2275 33 1 0.382 0.04905 0.2970 0.491
## 2285 32 1 0.370 0.04895 0.2855 0.480
## 2375 26 1 0.356 0.04909 0.2715 0.466
## 2437 21 1 0.339 0.04959 0.2544 0.451
## 2457 19 1 0.321 0.05009 0.2365 0.436
## 2476 17 1 0.302 0.05057 0.2176 0.419
## 2486 16 1 0.283 0.05082 0.1993 0.403
## 2541 14 1 0.263 0.05106 0.1798 0.385
## 2567 11 1 0.239 0.05171 0.1565 0.365
## 2583 10 1 0.215 0.05178 0.1343 0.345
## 2766 7 1 0.184 0.05272 0.1053 0.323
## 2823 6 1 0.154 0.05213 0.0791 0.299
##
## coxData$rltl_change_3G=3
## time n.risk n.event survival std.err lower 95% CI upper 95% CI
## 302 101 1 0.9901 0.00985 0.9710 1.000
## 495 100 1 0.9802 0.01386 0.9534 1.000
## 648 99 1 0.9703 0.01689 0.9377 1.000
## 658 98 1 0.9604 0.01941 0.9231 0.999
## 721 97 1 0.9505 0.02158 0.9091 0.994
## 735 96 1 0.9406 0.02352 0.8956 0.988
## 798 95 1 0.9307 0.02527 0.8825 0.982
## 803 94 1 0.9208 0.02687 0.8696 0.975
## 822 93 1 0.9109 0.02835 0.8570 0.968
## 830 92 1 0.9010 0.02972 0.8446 0.961
## 831 91 1 0.8911 0.03100 0.8324 0.954
## 833 90 1 0.8812 0.03220 0.8203 0.947
## 836 89 1 0.8713 0.03332 0.8084 0.939
## 840 88 1 0.8614 0.03438 0.7966 0.931
## 867 87 1 0.8515 0.03538 0.7849 0.924
## 895 86 1 0.8416 0.03633 0.7733 0.916
## 897 85 1 0.8317 0.03723 0.7618 0.908
## 919 84 1 0.8218 0.03808 0.7504 0.900
## 930 83 1 0.8119 0.03889 0.7391 0.892
## 932 82 1 0.8020 0.03965 0.7279 0.884
## 982 81 1 0.7921 0.04038 0.7168 0.875
## 1033 80 1 0.7822 0.04107 0.7057 0.867
## 1045 79 1 0.7723 0.04173 0.6947 0.859
## 1057 78 1 0.7624 0.04235 0.6837 0.850
## 1128 77 1 0.7525 0.04294 0.6728 0.842
## 1147 76 1 0.7426 0.04350 0.6620 0.833
## 1158 75 1 0.7327 0.04404 0.6513 0.824
## 1186 74 1 0.7228 0.04454 0.6405 0.816
## 1193 73 1 0.7129 0.04502 0.6299 0.807
## 1209 72 1 0.7030 0.04547 0.6193 0.798
## 1263 71 1 0.6931 0.04589 0.6087 0.789
## 1301 70 1 0.6832 0.04629 0.5982 0.780
## 1373 69 1 0.6733 0.04667 0.5877 0.771
## 1381 68 1 0.6634 0.04702 0.5773 0.762
## 1429 67 1 0.6535 0.04735 0.5669 0.753
## 1483 66 1 0.6436 0.04766 0.5566 0.744
## 1515 65 1 0.6337 0.04794 0.5463 0.735
## 1516 64 1 0.6238 0.04820 0.5361 0.726
## 1551 63 1 0.6139 0.04844 0.5259 0.717
## 1561 62 1 0.6040 0.04866 0.5157 0.707
## 1575 61 1 0.5941 0.04886 0.5056 0.698
## 1593 60 1 0.5842 0.04904 0.4955 0.689
## 1604 59 1 0.5743 0.04920 0.4855 0.679
## 1612 58 1 0.5644 0.04934 0.4755 0.670
## 1620 57 1 0.5545 0.04946 0.4655 0.660
## 1633 56 1 0.5446 0.04955 0.4556 0.651
## 1638 55 1 0.5347 0.04963 0.4457 0.641
## 1671 53 1 0.5246 0.04971 0.4356 0.632
## 1701 52 1 0.5145 0.04977 0.4256 0.622
## 1770 51 1 0.5044 0.04980 0.4156 0.612
## 1785 50 1 0.4943 0.04982 0.4057 0.602
## 1786 49 1 0.4842 0.04981 0.3958 0.592
## 1830 48 1 0.4741 0.04979 0.3859 0.582
## 1834 47 1 0.4640 0.04974 0.3761 0.573
## 1902 46 1 0.4540 0.04967 0.3663 0.563
## 1904 45 1 0.4439 0.04958 0.3566 0.552
## 1914 44 1 0.4338 0.04947 0.3469 0.542
## 1960 43 1 0.4237 0.04934 0.3372 0.532
## 1975 41 1 0.4134 0.04920 0.3273 0.522
## 1994 39 1 0.4028 0.04907 0.3172 0.511
## 2017 37 1 0.3919 0.04894 0.3068 0.501
## 2035 35 1 0.3807 0.04880 0.2961 0.489
## 2060 33 1 0.3691 0.04867 0.2851 0.478
## 2063 32 1 0.3576 0.04849 0.2741 0.466
## 2077 31 1 0.3461 0.04828 0.2633 0.455
## 2090 29 1 0.3341 0.04807 0.2520 0.443
## 2140 25 1 0.3208 0.04797 0.2393 0.430
## 2155 23 1 0.3068 0.04787 0.2260 0.417
## 2192 22 1 0.2929 0.04768 0.2129 0.403
## 2214 21 1 0.2789 0.04741 0.1999 0.389
## 2275 20 1 0.2650 0.04704 0.1871 0.375
## 2325 17 1 0.2494 0.04679 0.1727 0.360
## 2333 15 1 0.2328 0.04653 0.1573 0.344
## 2452 12 1 0.2134 0.04652 0.1392 0.327
## 2510 11 1 0.1940 0.04616 0.1217 0.309
## 2579 9 1 0.1724 0.04578 0.1025 0.290
## 2708 7 1 0.1478 0.04539 0.0810 0.270
## 2767 4 1 0.1108 0.04672 0.0485 0.253
## 2822 2 1 0.0554 0.04562 0.0110 0.278#### plot results for discrete scale
g4<-ggsurvplot(coxFit, data=coxData, pval = 0.003, ggtheme = theme_minimal(), risk.table.y.text.col = T, risk.table.y.text = FALSE, legend.labs=c("Tertile 1 (most RLTL attrition)", "Tertile 2", "Tertile 3 (least RLTL attrition)"))+theme_survminer(40)+ guides(colour = guide_legend(nrow = 3))
ggsurvplot(coxFit, data=coxData, pval = 0.003, ggtheme = theme_minimal(), risk.table.y.text.col = T, risk.table.y.text = FALSE, legend.labs=c("Tertile 1 (most RLTL attrition)", "Tertile 2", "Tertile 3 (least RLTL attrition)")) +theme_survminer(40)
#ggsave(file = "/Users/luise/Documents/Work/CattleTelomereLength/AgingCell/Images/RLTLchange.pdf", print(g4))ggplot(coxData, aes(x=as.factor(coxData$abs_rltl_change_3G), y= abs_rltl_change))+geom_boxplot() + xlab("Mean absolute RLTL change tertile") + ylab("Mean absolute RLTL change")+ theme_light(20)+ geom_jitter()## Warning: Use of `coxData$abs_rltl_change_3G` is discouraged. Use
## `abs_rltl_change_3G` instead.
## Warning: Use of `coxData$abs_rltl_change_3G` is discouraged. Use
## `abs_rltl_change_3G` instead.
ggplot(coxData, aes(x=as.factor(coxData$rltl_change_3G), y= rltl_change))+geom_boxplot()+ xlab("Mean RLTL change tertile") + ylab("Mean RLTL change") + theme_light(20)+ geom_jitter()## Warning: Use of `coxData$rltl_change_3G` is discouraged. Use `rltl_change_3G`
## instead.## Warning: Use of `coxData$rltl_change_3G` is discouraged. Use `rltl_change_3G`
## instead.
ggplot(coxData, aes(x=as.factor(coxData$average_resid_3G), y= average_resid))+geom_boxplot()+ xlab("Mean RLTL tertile") + ylab("Mean RLTL")+ theme_light(20)+ geom_jitter()## Warning: Use of `coxData$average_resid_3G` is discouraged. Use
## `average_resid_3G` instead.## Warning: Use of `coxData$average_resid_3G` is discouraged. Use
## `average_resid_3G` instead.
RLTL at the age of 1 year and productive lifespan
coxFit<-coxph(Surv(coxData$TimeMeasure, coxData$Event == 1)~coxData$av_rltl_1_res)
summary(coxFit)## Call:
## coxph(formula = Surv(coxData$TimeMeasure, coxData$Event == 1) ~
## coxData$av_rltl_1_res)
##
## n= 284, number of events= 221
## (21 observations deleted due to missingness)
##
## coef exp(coef) se(coef) z Pr(>|z|)
## coxData$av_rltl_1_res -0.9961 0.3693 0.4774 -2.087 0.0369 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## exp(coef) exp(-coef) lower .95 upper .95
## coxData$av_rltl_1_res 0.3693 2.708 0.1449 0.9414
##
## Concordance= 0.544 (se = 0.021 )
## Likelihood ratio test= 4.38 on 1 df, p=0.04
## Wald test = 4.35 on 1 df, p=0.04
## Score (logrank) test = 4.36 on 1 df, p=0.04All three measures of lifetime RLTL dynamics and productive lifespan
coxFit<-coxph(Surv(coxData$TimeMeasure, coxData$Event == 1)~coxData$rltl_change+ coxData$average_resid+ coxData$abs_rltl_change)
summary(coxFit)
Including average milk production
coxFit<-coxph(Surv(coxDataR$TimeMeasure, coxDataR$Event == 1)~coxDataR$average_resid+ as.numeric(paste(coxDataR$amp_305_lactation)) + coxDataR$totalTELnumber)
summary(coxFit)## Call:
## coxph(formula = Surv(coxDataR$TimeMeasure, coxDataR$Event ==
## 1) ~ coxDataR$average_resid + as.numeric(paste(coxDataR$amp_305_lactation)) +
## coxDataR$totalTELnumber)
##
## n= 253, number of events= 190
##
## coef exp(coef) se(coef)
## coxDataR$average_resid 6.860e-01 1.986e+00 6.655e-01
## as.numeric(paste(coxDataR$amp_305_lactation)) -6.233e-05 9.999e-01 2.878e-05
## coxDataR$totalTELnumber -6.125e-01 5.420e-01 9.017e-02
## z Pr(>|z|)
## coxDataR$average_resid 1.031 0.3026
## as.numeric(paste(coxDataR$amp_305_lactation)) -2.166 0.0303 *
## coxDataR$totalTELnumber -6.793 1.1e-11 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## exp(coef) exp(-coef) lower .95
## coxDataR$average_resid 1.9857 0.5036 0.5388
## as.numeric(paste(coxDataR$amp_305_lactation)) 0.9999 1.0001 0.9999
## coxDataR$totalTELnumber 0.5420 1.8451 0.4542
## upper .95
## coxDataR$average_resid 7.3177
## as.numeric(paste(coxDataR$amp_305_lactation)) 1.0000
## coxDataR$totalTELnumber 0.6467
##
## Concordance= 0.704 (se = 0.02 )
## Likelihood ratio test= 65.84 on 3 df, p=3e-14
## Wald test = 60.59 on 3 df, p=4e-13
## Score (logrank) test = 62.99 on 3 df, p=1e-13coxFit<-coxph(Surv(coxDataR$TimeMeasure, coxDataR$Event == 1)~coxDataR$abs_rltl_change+ as.numeric(paste(coxDataR$amp_305_lactation))+coxDataR$totalTELnumber)
summary(coxFit)## Call:
## coxph(formula = Surv(coxDataR$TimeMeasure, coxDataR$Event ==
## 1) ~ coxDataR$abs_rltl_change + as.numeric(paste(coxDataR$amp_305_lactation)) +
## coxDataR$totalTELnumber)
##
## n= 253, number of events= 190
##
## coef exp(coef) se(coef)
## coxDataR$abs_rltl_change 1.066e+00 2.905e+00 1.147e+00
## as.numeric(paste(coxDataR$amp_305_lactation)) -6.366e-05 9.999e-01 2.868e-05
## coxDataR$totalTELnumber -5.948e-01 5.517e-01 9.143e-02
## z Pr(>|z|)
## coxDataR$abs_rltl_change 0.930 0.3524
## as.numeric(paste(coxDataR$amp_305_lactation)) -2.219 0.0265 *
## coxDataR$totalTELnumber -6.506 7.73e-11 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## exp(coef) exp(-coef) lower .95
## coxDataR$abs_rltl_change 2.9046 0.3443 0.3070
## as.numeric(paste(coxDataR$amp_305_lactation)) 0.9999 1.0001 0.9999
## coxDataR$totalTELnumber 0.5517 1.8127 0.4612
## upper .95
## coxDataR$abs_rltl_change 27.4796
## as.numeric(paste(coxDataR$amp_305_lactation)) 1.0000
## coxDataR$totalTELnumber 0.6599
##
## Concordance= 0.703 (se = 0.02 )
## Likelihood ratio test= 65.62 on 3 df, p=4e-14
## Wald test = 60.45 on 3 df, p=5e-13
## Score (logrank) test = 63.14 on 3 df, p=1e-13coxFit<-coxph(Surv(coxDataR$TimeMeasure, coxDataR$Event == 1)~coxDataR$rltl_change+ as.numeric(paste(coxDataR$amp_305_lactation))+coxDataR$totalTELnumber)
summary(coxFit)## Call:
## coxph(formula = Surv(coxDataR$TimeMeasure, coxDataR$Event ==
## 1) ~ coxDataR$rltl_change + as.numeric(paste(coxDataR$amp_305_lactation)) +
## coxDataR$totalTELnumber)
##
## n= 253, number of events= 190
##
## coef exp(coef) se(coef)
## coxDataR$rltl_change -3.301e+00 3.685e-02 1.214e+00
## as.numeric(paste(coxDataR$amp_305_lactation)) -5.700e-05 9.999e-01 2.859e-05
## coxDataR$totalTELnumber -5.761e-01 5.621e-01 9.059e-02
## z Pr(>|z|)
## coxDataR$rltl_change -2.719 0.00655 **
## as.numeric(paste(coxDataR$amp_305_lactation)) -1.994 0.04617 *
## coxDataR$totalTELnumber -6.359 2.03e-10 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## exp(coef) exp(-coef) lower .95
## coxDataR$rltl_change 0.03685 27.140 0.003411
## as.numeric(paste(coxDataR$amp_305_lactation)) 0.99994 1.000 0.999887
## coxDataR$totalTELnumber 0.56211 1.779 0.470665
## upper .95
## coxDataR$rltl_change 0.3980
## as.numeric(paste(coxDataR$amp_305_lactation)) 1.0000
## coxDataR$totalTELnumber 0.6713
##
## Concordance= 0.706 (se = 0.02 )
## Likelihood ratio test= 71.88 on 3 df, p=2e-15
## Wald test = 70.78 on 3 df, p=3e-15
## Score (logrank) test = 72.11 on 3 df, p=2e-15coxFit<-coxph(Surv(coxDataR$TimeMeasure, coxDataR$Event == 1)~coxDataR$av_rltl_1_res + as.numeric(paste(coxDataR$amp_305_lactation))+ coxDataR$totalTELnumber)
summary(coxFit)## Call:
## coxph(formula = Surv(coxDataR$TimeMeasure, coxDataR$Event ==
## 1) ~ coxDataR$av_rltl_1_res + as.numeric(paste(coxDataR$amp_305_lactation)) +
## coxDataR$totalTELnumber)
##
## n= 239, number of events= 177
## (14 observations deleted due to missingness)
##
## coef exp(coef) se(coef)
## coxDataR$av_rltl_1_res -4.681e-02 9.543e-01 5.445e-01
## as.numeric(paste(coxDataR$amp_305_lactation)) -5.465e-05 9.999e-01 2.984e-05
## coxDataR$totalTELnumber -5.417e-01 5.818e-01 9.348e-02
## z Pr(>|z|)
## coxDataR$av_rltl_1_res -0.086 0.9315
## as.numeric(paste(coxDataR$amp_305_lactation)) -1.831 0.0671 .
## coxDataR$totalTELnumber -5.794 6.86e-09 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## exp(coef) exp(-coef) lower .95
## coxDataR$av_rltl_1_res 0.9543 1.048 0.3282
## as.numeric(paste(coxDataR$amp_305_lactation)) 0.9999 1.000 0.9999
## coxDataR$totalTELnumber 0.5818 1.719 0.4844
## upper .95
## coxDataR$av_rltl_1_res 2.7745
## as.numeric(paste(coxDataR$amp_305_lactation)) 1.0000
## coxDataR$totalTELnumber 0.6988
##
## Concordance= 0.686 (se = 0.021 )
## Likelihood ratio test= 48.63 on 3 df, p=2e-10
## Wald test = 44.57 on 3 df, p=1e-09
## Score (logrank) test = 46.14 on 3 df, p=5e-10#total TEL sample number was not significant and therefore will be eliminated:I THINK I SHOULD USE THE MODELS BELOW
coxFit<-coxph(Surv(coxDataR$TimeMeasure, coxDataR$Event == 1)~coxDataR$average_resid+ as.numeric(paste(coxDataR$amp_305_lactation)))
summary(coxFit)## Call:
## coxph(formula = Surv(coxDataR$TimeMeasure, coxDataR$Event ==
## 1) ~ coxDataR$average_resid + as.numeric(paste(coxDataR$amp_305_lactation)))
##
## n= 253, number of events= 190
##
## coef exp(coef) se(coef)
## coxDataR$average_resid 2.315e-01 1.260e+00 6.708e-01
## as.numeric(paste(coxDataR$amp_305_lactation)) -9.933e-05 9.999e-01 2.822e-05
## z Pr(>|z|)
## coxDataR$average_resid 0.345 0.730011
## as.numeric(paste(coxDataR$amp_305_lactation)) -3.520 0.000432 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## exp(coef) exp(-coef) lower .95
## coxDataR$average_resid 1.2605 0.7934 0.3385
## as.numeric(paste(coxDataR$amp_305_lactation)) 0.9999 1.0001 0.9998
## upper .95
## coxDataR$average_resid 4.693
## as.numeric(paste(coxDataR$amp_305_lactation)) 1.000
##
## Concordance= 0.598 (se = 0.023 )
## Likelihood ratio test= 11.93 on 2 df, p=0.003
## Wald test = 12.43 on 2 df, p=0.002
## Score (logrank) test = 12.44 on 2 df, p=0.002coxFit<-coxph(Surv(coxDataR$TimeMeasure, coxDataR$Event == 1)~coxDataR$abs_rltl_change+ as.numeric(paste(coxDataR$amp_305_lactation)))
summary(coxFit)## Call:
## coxph(formula = Surv(coxDataR$TimeMeasure, coxDataR$Event ==
## 1) ~ coxDataR$abs_rltl_change + as.numeric(paste(coxDataR$amp_305_lactation)))
##
## n= 253, number of events= 190
##
## coef exp(coef) se(coef)
## coxDataR$abs_rltl_change 2.411e+00 1.115e+01 1.146e+00
## as.numeric(paste(coxDataR$amp_305_lactation)) -1.018e-04 9.999e-01 2.818e-05
## z Pr(>|z|)
## coxDataR$abs_rltl_change 2.104 0.035391 *
## as.numeric(paste(coxDataR$amp_305_lactation)) -3.611 0.000305 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## exp(coef) exp(-coef) lower .95
## coxDataR$abs_rltl_change 11.1466 0.08971 1.1793
## as.numeric(paste(coxDataR$amp_305_lactation)) 0.9999 1.00010 0.9998
## upper .95
## coxDataR$abs_rltl_change 105.4
## as.numeric(paste(coxDataR$amp_305_lactation)) 1.0
##
## Concordance= 0.598 (se = 0.023 )
## Likelihood ratio test= 15.99 on 2 df, p=3e-04
## Wald test = 17.06 on 2 df, p=2e-04
## Score (logrank) test = 17.08 on 2 df, p=2e-04coxFit<-coxph(Surv(coxDataR$TimeMeasure, coxDataR$Event == 1)~coxDataR$rltl_change+ as.numeric(paste(coxDataR$amp_305_lactation)))
summary(coxFit)## Call:
## coxph(formula = Surv(coxDataR$TimeMeasure, coxDataR$Event ==
## 1) ~ coxDataR$rltl_change + as.numeric(paste(coxDataR$amp_305_lactation)))
##
## n= 253, number of events= 190
##
## coef exp(coef) se(coef)
## coxDataR$rltl_change -5.035e+00 6.506e-03 1.315e+00
## as.numeric(paste(coxDataR$amp_305_lactation)) -9.057e-05 9.999e-01 2.819e-05
## z Pr(>|z|)
## coxDataR$rltl_change -3.829 0.000129 ***
## as.numeric(paste(coxDataR$amp_305_lactation)) -3.212 0.001317 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## exp(coef) exp(-coef) lower .95
## coxDataR$rltl_change 0.006506 153.7 0.0004944
## as.numeric(paste(coxDataR$amp_305_lactation)) 0.999909 1.0 0.9998542
## upper .95
## coxDataR$rltl_change 0.08562
## as.numeric(paste(coxDataR$amp_305_lactation)) 0.99996
##
## Concordance= 0.601 (se = 0.022 )
## Likelihood ratio test= 25.19 on 2 df, p=3e-06
## Wald test = 28.18 on 2 df, p=8e-07
## Score (logrank) test = 26.7 on 2 df, p=2e-06test RLTL at the age of 1 year in the same model as mean RLTL change
coxFit<-coxph(Surv(coxDataR$TimeMeasure, coxDataR$Event == 1)~coxDataR$rltl_change+ coxDataR$av_rltl_1_res+ as.numeric(paste(coxDataR$amp_305_lactation)))
summary(coxFit)## Call:
## coxph(formula = Surv(coxDataR$TimeMeasure, coxDataR$Event ==
## 1) ~ coxDataR$rltl_change + coxDataR$av_rltl_1_res + as.numeric(paste(coxDataR$amp_305_lactation)))
##
## n= 239, number of events= 177
## (14 observations deleted due to missingness)
##
## coef exp(coef) se(coef)
## coxDataR$rltl_change -5.030e+00 6.537e-03 1.379e+00
## coxDataR$av_rltl_1_res -6.552e-01 5.193e-01 5.263e-01
## as.numeric(paste(coxDataR$amp_305_lactation)) -7.896e-05 9.999e-01 2.928e-05
## z Pr(>|z|)
## coxDataR$rltl_change -3.648 0.000264 ***
## coxDataR$av_rltl_1_res -1.245 0.213134
## as.numeric(paste(coxDataR$amp_305_lactation)) -2.696 0.007008 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## exp(coef) exp(-coef) lower .95
## coxDataR$rltl_change 0.006537 152.972 0.0004382
## coxDataR$av_rltl_1_res 0.519327 1.926 0.1851255
## as.numeric(paste(coxDataR$amp_305_lactation)) 0.999921 1.000 0.9998637
## upper .95
## coxDataR$rltl_change 0.09752
## coxDataR$av_rltl_1_res 1.45685
## as.numeric(paste(coxDataR$amp_305_lactation)) 0.99998
##
## Concordance= 0.595 (se = 0.023 )
## Likelihood ratio test= 22.47 on 3 df, p=5e-05
## Wald test = 25.18 on 3 df, p=1e-05
## Score (logrank) test = 23.85 on 3 df, p=3e-05All three measures of lifetime RLTL dynamics in the same model including milk productivity measure:
coxFit<-coxph(Surv(coxDataR$TimeMeasure, coxDataR$Event == 1)~coxDataR$rltl_change+coxDataR$average_resid+coxDataR$abs_rltl_change+ as.numeric(paste(coxDataR$amp_305_lactation)))
summary(coxFit)## Call:
## coxph(formula = Surv(coxDataR$TimeMeasure, coxDataR$Event ==
## 1) ~ coxDataR$rltl_change + coxDataR$average_resid + coxDataR$abs_rltl_change +
## as.numeric(paste(coxDataR$amp_305_lactation)))
##
## n= 253, number of events= 190
##
## coef exp(coef) se(coef)
## coxDataR$rltl_change -4.525e+00 1.083e-02 1.387e+00
## coxDataR$average_resid -6.968e-02 9.327e-01 6.840e-01
## coxDataR$abs_rltl_change 1.237e+00 3.447e+00 1.224e+00
## as.numeric(paste(coxDataR$amp_305_lactation)) -9.162e-05 9.999e-01 2.824e-05
## z Pr(>|z|)
## coxDataR$rltl_change -3.262 0.00111 **
## coxDataR$average_resid -0.102 0.91885
## coxDataR$abs_rltl_change 1.011 0.31198
## as.numeric(paste(coxDataR$amp_305_lactation)) -3.244 0.00118 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## exp(coef) exp(-coef) lower .95
## coxDataR$rltl_change 0.01083 92.3027 0.0007142
## coxDataR$average_resid 0.93269 1.0722 0.2440869
## coxDataR$abs_rltl_change 3.44666 0.2901 0.3130755
## as.numeric(paste(coxDataR$amp_305_lactation)) 0.99991 1.0001 0.9998530
## upper .95
## coxDataR$rltl_change 0.1643
## coxDataR$average_resid 3.5639
## coxDataR$abs_rltl_change 37.9445
## as.numeric(paste(coxDataR$amp_305_lactation)) 1.0000
##
## Concordance= 0.6 (se = 0.022 )
## Likelihood ratio test= 26.21 on 4 df, p=3e-05
## Wald test = 30.41 on 4 df, p=4e-06
## Score (logrank) test = 28.88 on 4 df, p=8e-06Repeat Cox proportional hazard model without without the fist measurement.
red_data_l_new<-read.csv("/Users/luise/Documents/Work/CattleTelomereLength/LongTelomereChangePaper/git/TelomereChangeInDairyCattle/excluded_less_3_month.csv")uniq_red_data_l_new<-red_data_l_new[!duplicated(red_data_l_new$recoded_id),]Histogram of samples per animal after removing early measurements
hist1b<-ggplot(uniq_red_data_l_new, aes(x = uniq_red_data_l_new$SamplesPerAnimal)) +
geom_histogram(colour= "black", fill = mycoloursP[15], binwidth = 1) + xlab("Number of samples per animal") +
geom_vline(xintercept= mean(as.numeric(paste(uniq_red_data_l_new$SamplesPerAnimal)), na.rm=TRUE), colour="red")+
annotate("text", x = 0.7, y = 100, label = "(A)", fontface =2, size = 10) + theme_classic(20)+
scale_x_continuous(breaks=c(1:8))## Warning in mean(as.numeric(paste(uniq_red_data_l_new$SamplesPerAnimal)), : NAs
## introduced by coercionhist1b## Warning: Use of `uniq_red_data_l_new$SamplesPerAnimal` is discouraged. Use
## `SamplesPerAnimal` instead.## Warning: Removed 27 rows containing non-finite values (stat_bin).
For dataframe that excludes the first sample
hist2b<-ggplot(red_data_l_new, aes(x = red_data_l_new$residual_rltl_t)) +
geom_histogram(colour= "black", fill = mycoloursP[16]) + xlab("Adjusted RLTL") +
geom_vline(xintercept= mean(as.numeric(paste(red_data_l_new$residual_rltl_t))), colour="red") +
annotate("text", x = -0.45, y = 120, label = "(B)", fontface =2, size = 10)+
theme_classic(20)+
stat_function(
fun = function(x, mean, sd, n){
n/27* dnorm(x = x, mean = mean, sd = sd)
},
args = with(red_data_l_new, c(mean = mean(red_data_l_new$residual_rltl_t), sd = sd(red_data_l_new$residual_rltl_t), n
= length(red_data_l_new$residual_rltl_t)))
)
hist2b## Warning: Use of `red_data_l_new$residual_rltl_t` is discouraged. Use
## `residual_rltl_t` instead.
## Warning: Use of `red_data_l_new$residual_rltl_t` is discouraged. Use
## `residual_rltl_t` instead.## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
#red_data_l_new_ch<-subset(red_data_l_new, red_data_l_new$SamplesPerAnimal != 2)
diff<- as.numeric(paste(red_data_l_new$diff_rltl_res))
hist3b<-ggplot(red_data_l_new, aes(x = as.numeric(paste(red_data_l_new$diff_rltl_res)))) +
geom_histogram(colour= "black", fill = mycoloursP[17]) + xlab("Adjusted RLTL change") +
geom_vline(xintercept= mean(as.numeric(paste(red_data_l_new$diff_rltl_res))), colour="red") +
annotate("text", x = -0.75, y = 95, label = "(C)", fontface =2, size= 10)+
theme_classic(20)+
stat_function(
fun = function(x, mean, sd, n){
n/23* dnorm(x = x, mean = mean, sd = sd)
},
args = with(red_data_l_new, c(mean = mean(as.numeric(paste(red_data_l_new$diff_rltl_res))), sd = sd(as.numeric(paste(red_data_l_new$diff_rltl_res))), n
= length(as.numeric(paste(red_data_l_new$diff_rltl_res))))))
hist3b## Warning: Use of `red_data_l_new$diff_rltl_res` is discouraged. Use
## `diff_rltl_res` instead.
## Warning: Use of `red_data_l_new$diff_rltl_res` is discouraged. Use
## `diff_rltl_res` instead.## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
The plot above looks so much like hist1 that I think I may have made a mistake. when I look for summary output of diff_rltl_res in both the complete dataset (data) as well as the subsetted dataset (red_data_l_new) summary statistics are really very similar, but not identical. By removing single samples, I should have removed samples without a first change measurement. I have to think about this for a little bit.
summary(as.numeric(paste(red_data_l_new$diff_rltl_res)))## Min. 1st Qu. Median Mean 3rd Qu. Max.
## -0.72730 -0.15316 -0.03555 -0.03257 0.08741 0.48019summary(as.numeric(paste(data$diff_rltl_res)))## Warning in summary(as.numeric(paste(data$diff_rltl_res))): NAs introduced by
## coercion## Min. 1st Qu. Median Mean 3rd Qu. Max. NA's
## -0.72730 -0.15310 -0.03538 -0.03243 0.08796 0.48019 305u_red_data_l_new<-red_data_l_new[!duplicated(red_data_l_new$recoded_id),]
coxDataL<-data.frame(recoded_id=u_red_data_l_new$recoded_id, dob=u_red_data_l_new$dob, herd_life=u_red_data_l_new$herd_life, cull_date=u_red_data_l_new$cull_date, average_resid=u_red_data_l_new$average_resid, rltl_change=u_red_data_l_new$rltl_change, abs_rltl_change=u_red_data_l_new$abs_rltl_change, amp_305_lactation= u_red_data_l_new$amp_305_lactation)
coxDataL$CutoffDate<-"01/01/2017"
coxDataL$Event<-ifelse(coxDataL$cull_date=="NULL", 0, 1)
coxDataL$TimeMeasure<-ifelse(coxDataL$Event==0, as.Date(coxDataL$CutoffDate, "%d/%m/%Y")- as.Date(coxDataL$dob, "%d/%m/%Y"), as.numeric(paste(coxDataL$herd_life)))## Warning in ifelse(coxDataL$Event == 0, as.Date(coxDataL$CutoffDate, "%d/%m/%Y")
## - : NAs introduced by coercioncoxFit<-coxph(Surv(coxDataL$TimeMeasure, coxDataL$Event == 1)~coxDataL$average_resid)
summary(coxFit)## Call:
## coxph(formula = Surv(coxDataL$TimeMeasure, coxDataL$Event ==
## 1) ~ coxDataL$average_resid)
##
## n= 304, number of events= 240
##
## coef exp(coef) se(coef) z Pr(>|z|)
## coxDataL$average_resid -0.5647 0.5686 0.5667 -0.996 0.319
##
## exp(coef) exp(-coef) lower .95 upper .95
## coxDataL$average_resid 0.5686 1.759 0.1873 1.726
##
## Concordance= 0.524 (se = 0.021 )
## Likelihood ratio test= 1 on 1 df, p=0.3
## Wald test = 0.99 on 1 df, p=0.3
## Score (logrank) test = 0.99 on 1 df, p=0.3coxFit<-coxph(Surv(coxDataL$TimeMeasure, coxDataL$Event == 1)~coxDataL$abs_rltl_change)
summary(coxFit)## Call:
## coxph(formula = Surv(coxDataL$TimeMeasure, coxDataL$Event ==
## 1) ~ coxDataL$abs_rltl_change)
##
## n= 304, number of events= 240
##
## coef exp(coef) se(coef) z Pr(>|z|)
## coxDataL$abs_rltl_change 2.7934 16.3372 0.9814 2.846 0.00442 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## exp(coef) exp(-coef) lower .95 upper .95
## coxDataL$abs_rltl_change 16.34 0.06121 2.387 111.8
##
## Concordance= 0.537 (se = 0.023 )
## Likelihood ratio test= 7.49 on 1 df, p=0.006
## Wald test = 8.1 on 1 df, p=0.004
## Score (logrank) test = 8.05 on 1 df, p=0.005coxFit<-coxph(Surv(coxDataL$TimeMeasure, coxDataL$Event == 1)~coxDataL$rltl_change)
summary(coxFit)## Call:
## coxph(formula = Surv(coxDataL$TimeMeasure, coxDataL$Event ==
## 1) ~ coxDataL$rltl_change)
##
## n= 304, number of events= 240
##
## coef exp(coef) se(coef) z Pr(>|z|)
## coxDataL$rltl_change -5.065332 0.006312 0.869793 -5.824 5.76e-09 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## exp(coef) exp(-coef) lower .95 upper .95
## coxDataL$rltl_change 0.006312 158.4 0.001148 0.03472
##
## Concordance= 0.572 (se = 0.023 )
## Likelihood ratio test= 26.01 on 1 df, p=3e-07
## Wald test = 33.91 on 1 df, p=6e-09
## Score (logrank) test = 30.05 on 1 df, p=4e-08#now without first measurement but including milk productivity
#now I will include average milk production, TEL sample number was tested, but removed from the models because it was not statistically associated with productive lifespan
coxDataRL<-subset(coxDataL, coxDataL$amp_305_lactation != "NULL")
#coxDataRL<-merge(coxDataRL, maxSN, by = "recoded_id")
coxFit<-coxph(Surv(coxDataRL$TimeMeasure, coxDataRL$Event == 1)~coxDataRL$average_resid+ as.numeric(paste(coxDataRL$amp_305_lactation)))
summary(coxFit)## Call:
## coxph(formula = Surv(coxDataRL$TimeMeasure, coxDataRL$Event ==
## 1) ~ coxDataRL$average_resid + as.numeric(paste(coxDataRL$amp_305_lactation)))
##
## n= 253, number of events= 190
##
## coef exp(coef) se(coef)
## coxDataRL$average_resid -4.847e-01 6.159e-01 6.494e-01
## as.numeric(paste(coxDataRL$amp_305_lactation)) -9.853e-05 9.999e-01 2.812e-05
## z Pr(>|z|)
## coxDataRL$average_resid -0.746 0.455445
## as.numeric(paste(coxDataRL$amp_305_lactation)) -3.504 0.000458 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## exp(coef) exp(-coef) lower .95
## coxDataRL$average_resid 0.6159 1.624 0.1725
## as.numeric(paste(coxDataRL$amp_305_lactation)) 0.9999 1.000 0.9998
## upper .95
## coxDataRL$average_resid 2.199
## as.numeric(paste(coxDataRL$amp_305_lactation)) 1.000
##
## Concordance= 0.599 (se = 0.023 )
## Likelihood ratio test= 12.37 on 2 df, p=0.002
## Wald test = 12.93 on 2 df, p=0.002
## Score (logrank) test = 12.92 on 2 df, p=0.002coxFit<-coxph(Surv(coxDataRL$TimeMeasure, coxDataRL$Event == 1)~coxDataRL$abs_rltl_change+ as.numeric(paste(coxDataRL$amp_305_lactation)))
summary(coxFit)## Call:
## coxph(formula = Surv(coxDataRL$TimeMeasure, coxDataRL$Event ==
## 1) ~ coxDataRL$abs_rltl_change + as.numeric(paste(coxDataRL$amp_305_lactation)))
##
## n= 253, number of events= 190
##
## coef exp(coef) se(coef)
## coxDataRL$abs_rltl_change 2.403e+00 1.105e+01 1.147e+00
## as.numeric(paste(coxDataRL$amp_305_lactation)) -1.018e-04 9.999e-01 2.818e-05
## z Pr(>|z|)
## coxDataRL$abs_rltl_change 2.095 0.036134 *
## as.numeric(paste(coxDataRL$amp_305_lactation)) -3.610 0.000306 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## exp(coef) exp(-coef) lower .95
## coxDataRL$abs_rltl_change 11.0524 0.09048 1.1680
## as.numeric(paste(coxDataRL$amp_305_lactation)) 0.9999 1.00010 0.9998
## upper .95
## coxDataRL$abs_rltl_change 104.6
## as.numeric(paste(coxDataRL$amp_305_lactation)) 1.0
##
## Concordance= 0.598 (se = 0.023 )
## Likelihood ratio test= 15.95 on 2 df, p=3e-04
## Wald test = 17.02 on 2 df, p=2e-04
## Score (logrank) test = 17.04 on 2 df, p=2e-04coxFit<-coxph(Surv(coxDataRL$TimeMeasure, coxDataRL$Event == 1)~coxDataRL$rltl_change + as.numeric(paste(coxDataRL$amp_305_lactation)))
summary(coxFit)## Call:
## coxph(formula = Surv(coxDataRL$TimeMeasure, coxDataRL$Event ==
## 1) ~ coxDataRL$rltl_change + as.numeric(paste(coxDataRL$amp_305_lactation)))
##
## n= 253, number of events= 190
##
## coef exp(coef) se(coef)
## coxDataRL$rltl_change -5.056e+00 6.370e-03 1.315e+00
## as.numeric(paste(coxDataRL$amp_305_lactation)) -9.096e-05 9.999e-01 2.818e-05
## z Pr(>|z|)
## coxDataRL$rltl_change -3.844 0.000121 ***
## as.numeric(paste(coxDataRL$amp_305_lactation)) -3.228 0.001247 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## exp(coef) exp(-coef) lower .95
## coxDataRL$rltl_change 0.00637 157 0.0004836
## as.numeric(paste(coxDataRL$amp_305_lactation)) 0.99991 1 0.9998538
## upper .95
## coxDataRL$rltl_change 0.0839
## as.numeric(paste(coxDataRL$amp_305_lactation)) 1.0000
##
## Concordance= 0.601 (se = 0.022 )
## Likelihood ratio test= 25.3 on 2 df, p=3e-06
## Wald test = 28.32 on 2 df, p=7e-07
## Score (logrank) test = 26.79 on 2 df, p=2e-06coxFit<-coxph(Surv(coxDataRL$TimeMeasure, coxDataRL$Event == 1)~coxDataRL$rltl_change + coxDataRL$abs_rltl_change+ coxDataRL$average_resid+ as.numeric(paste(coxDataRL$amp_305_lactation)))
summary(coxFit)
See if relationship remains statistically significant when only animals with more than 2 samples are included in the analysis
coxData2<-data.frame(recoded_id=u_red_more_2$recoded_id, dob = u_red_more_2$dob, herd_life = u_red_more_2$herd_life, cull_date = u_red_more_2$cull_date, average_resid_M2 = u_red_more_2$average_resid_M2, rltl_change_M2 = u_red_more_2$rltl_change_M2, abs_rltl_change_M2 = u_red_more_2$abs_rltl_change_M2, amp_305_lactation = u_red_more_2$amp_305_lactation, av_rltl_1_res = u_red_more_2$av_rltl_1_res)
coxData2$CutoffDate<-"01/01/2017"
coxData2$Event<-ifelse(coxData2$cull_date=="NULL", 0, 1)
coxData2$TimeMeasure<-ifelse(coxData2$Event==0, as.Date(coxData2$CutoffDate, "%d/%m/%Y")- as.Date(coxData2$dob, "%d/%m/%Y"), as.numeric(paste(coxData2$herd_life)))## Warning in ifelse(coxData2$Event == 0, as.Date(coxData2$CutoffDate, "%d/%m/%Y")
## - : NAs introduced by coercioncoxData2$average_resid_3G<-ntile(coxData2$average_resid_M2,3)
coxData2$rltl_change_3G<-ntile(coxData2$rltl_change_M2,3)
coxData2$abs_rltl_change_3G<-ntile(coxData2$abs_rltl_change_M2,3)
# models:
##Average lifetime telomere length
### condinuous variable:
coxFit<-coxph(Surv(coxData2$TimeMeasure, coxData2$Event == 1)~coxData2$average_resid_M2)
summary(coxFit)## Call:
## coxph(formula = Surv(coxData2$TimeMeasure, coxData2$Event ==
## 1) ~ coxData2$average_resid_M2)
##
## n= 277, number of events= 213
##
## coef exp(coef) se(coef) z Pr(>|z|)
## coxData2$average_resid_M2 -0.5990 0.5494 0.6041 -0.992 0.321
##
## exp(coef) exp(-coef) lower .95 upper .95
## coxData2$average_resid_M2 0.5494 1.82 0.1682 1.795
##
## Concordance= 0.524 (se = 0.022 )
## Likelihood ratio test= 0.99 on 1 df, p=0.3
## Wald test = 0.98 on 1 df, p=0.3
## Score (logrank) test = 0.98 on 1 df, p=0.3### discrete scale
coxFit<-coxph(Surv(coxData2$TimeMeasure, coxData2$Event == 1)~coxData2$average_resid_3G)
summary(coxFit)## Call:
## coxph(formula = Surv(coxData2$TimeMeasure, coxData2$Event ==
## 1) ~ coxData2$average_resid_3G)
##
## n= 277, number of events= 213
##
## coef exp(coef) se(coef) z Pr(>|z|)
## coxData2$average_resid_3G -0.10350 0.90168 0.08528 -1.214 0.225
##
## exp(coef) exp(-coef) lower .95 upper .95
## coxData2$average_resid_3G 0.9017 1.109 0.7629 1.066
##
## Concordance= 0.526 (se = 0.021 )
## Likelihood ratio test= 1.47 on 1 df, p=0.2
## Wald test = 1.47 on 1 df, p=0.2
## Score (logrank) test = 1.48 on 1 df, p=0.2coxFit<-survfit(Surv(coxData2$TimeMeasure, coxData2$Event == 1)~coxData2$average_resid_3G)#### plot results for discrete scale
survp<- ggsurvplot(coxFit, data=coxData2, pval = TRUE, ggtheme = theme_minimal(), risk.table.y.text.col = T, risk.table.y.text = FALSE, legend.labs=c("Tertile 1 (shortest mean RLTL)", "Tertile 2", "Tertile 3 (longest mean RLTL)"))+theme_survminer(40)
survp
#ggsave(file = "/Users/luise/Documents/Work/CattleTelomereLength/AgingCell/Images/AverageRLTL.pdf", print(survp))Mean RLTL change
coxFit<-coxph(Surv(coxData2$TimeMeasure, coxData2$Event == 1)~coxData2$rltl_change_M2)
summary(coxFit)## Call:
## coxph(formula = Surv(coxData2$TimeMeasure, coxData2$Event ==
## 1) ~ coxData2$rltl_change_M2)
##
## n= 277, number of events= 213
##
## coef exp(coef) se(coef) z Pr(>|z|)
## coxData2$rltl_change_M2 -3.46163 0.03138 1.36042 -2.545 0.0109 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## exp(coef) exp(-coef) lower .95 upper .95
## coxData2$rltl_change_M2 0.03138 31.87 0.002181 0.4515
##
## Concordance= 0.544 (se = 0.023 )
## Likelihood ratio test= 6.42 on 1 df, p=0.01
## Wald test = 6.47 on 1 df, p=0.01
## Score (logrank) test = 6.46 on 1 df, p=0.01### discrete scale
coxFit<-coxph(Surv(coxData2$TimeMeasure, coxData2$Event == 1)~coxData2$rltl_change_3G)
summary(coxFit)## Call:
## coxph(formula = Surv(coxData2$TimeMeasure, coxData2$Event ==
## 1) ~ coxData2$rltl_change_3G)
##
## n= 277, number of events= 213
##
## coef exp(coef) se(coef) z Pr(>|z|)
## coxData2$rltl_change_3G -0.14184 0.86776 0.08938 -1.587 0.113
##
## exp(coef) exp(-coef) lower .95 upper .95
## coxData2$rltl_change_3G 0.8678 1.152 0.7283 1.034
##
## Concordance= 0.534 (se = 0.021 )
## Likelihood ratio test= 2.52 on 1 df, p=0.1
## Wald test = 2.52 on 1 df, p=0.1
## Score (logrank) test = 2.52 on 1 df, p=0.1coxFit<-survfit(Surv(coxData2$TimeMeasure, coxData2$Event == 1)~coxData2$rltl_change_3G)#### plot results for discrete scale
survp<- ggsurvplot(coxFit, data=coxData2, pval = TRUE, ggtheme = theme_minimal(), risk.table.y.text.col = T, risk.table.y.text = FALSE, legend.labs=c("Tertile 1 (most RLTL shortening)", "Tertile 2", "Tertile 3 (least RLTL shortening)"))+theme_survminer(40)
survp
mean absolute RLTL change
coxFit<-coxph(Surv(coxData2$TimeMeasure, coxData2$Event == 1)~coxData2$abs_rltl_change_M2)
summary(coxFit)## Call:
## coxph(formula = Surv(coxData2$TimeMeasure, coxData2$Event ==
## 1) ~ coxData2$abs_rltl_change_M2)
##
## n= 277, number of events= 213
##
## coef exp(coef) se(coef) z Pr(>|z|)
## coxData2$abs_rltl_change_M2 1.522 4.579 1.159 1.312 0.189
##
## exp(coef) exp(-coef) lower .95 upper .95
## coxData2$abs_rltl_change_M2 4.579 0.2184 0.4719 44.44
##
## Concordance= 0.528 (se = 0.024 )
## Likelihood ratio test= 1.67 on 1 df, p=0.2
## Wald test = 1.72 on 1 df, p=0.2
## Score (logrank) test = 1.72 on 1 df, p=0.2### discrete scale
coxFit<-coxph(Surv(coxData2$TimeMeasure, coxData2$Event == 1)~coxData2$abs_rltl_change_3G)
summary(coxFit)## Call:
## coxph(formula = Surv(coxData2$TimeMeasure, coxData2$Event ==
## 1) ~ coxData2$abs_rltl_change_3G)
##
## n= 277, number of events= 213
##
## coef exp(coef) se(coef) z Pr(>|z|)
## coxData2$abs_rltl_change_3G 0.16969 1.18494 0.08624 1.968 0.0491 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## exp(coef) exp(-coef) lower .95 upper .95
## coxData2$abs_rltl_change_3G 1.185 0.8439 1.001 1.403
##
## Concordance= 0.537 (se = 0.021 )
## Likelihood ratio test= 3.88 on 1 df, p=0.05
## Wald test = 3.87 on 1 df, p=0.05
## Score (logrank) test = 3.89 on 1 df, p=0.05coxFit<-survfit(Surv(coxData2$TimeMeasure, coxData2$Event == 1)~coxData2$abs_rltl_change_3G)#### plot results for discrete scale
survp<- ggsurvplot(coxFit, data=coxData2, pval = TRUE, ggtheme = theme_minimal(), risk.table.y.text.col = T, risk.table.y.text = FALSE, legend.labs=c("Tertile 1 (least absolute RLTL change)", "Tertile 2", "Tertile 3 (most absolute RLTL change)"))+theme_survminer(40)
survp
Some more figures:
variable<-as.numeric(paste(data$sample_interval))## Warning: NAs introduced by coercionmeanVar<-mean(variable, na.rm= TRUE)
sdVar<-sd(variable, na.rm=TRUE)
ggplot(data, aes(x = variable)) +
geom_histogram(colour= "black", fill = mycoloursP[7]) + xlab("Sample interval in days") + geom_vline(xintercept= meanVar, colour="red")+
theme_classic(20)## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.## Warning: Removed 305 rows containing non-finite values (stat_bin).
uData<-data[!duplicated(data$recoded_id),]
variable<-as.factor(paste(uData$birth_year))
ggplot(uData, aes(x = variable)) +
geom_bar(colour= "black", fill = mycoloursP[8]) + xlab("Birth year")+ theme_classic(20)
#Average milk production
variable<-as.numeric(paste(uData$amp_305_lactation))## Warning: NAs introduced by coercionmeanVar<-mean(variable, na.rm= TRUE)
sdVar<-sd(variable, na.rm=TRUE)
ggplot(uData, aes(x = variable)) +
geom_histogram(colour= "black", fill = mycoloursP[9]) + xlab("Average milk production") + geom_vline(xintercept= meanVar, colour="red")+
theme_classic(20)## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.## Warning: Removed 52 rows containing non-finite values (stat_bin).
# cull reasons
variable<-as.factor(paste(uDeadData$cull_reason))
ggplot(uDeadData, aes(x = variable)) +
geom_bar(colour= "black", fill = mycoloursP[10]) + xlab("Cull reason") +
theme_classic(10)+
theme(axis.text.x = element_text(angle = 90, hjust = 1))
# fertility events
newdata <- subset(uData, uData$health_event_fertil_all != "NULL")
variable<-as.integer(paste(newdata$health_event_fertil_all))
ggplot(newdata, aes(x = variable)) +
geom_bar(colour= "black", fill = mycoloursP[11]) + xlab("Fertility events")+
theme_classic(20)+
scale_x_continuous(breaks=c(0:12))
#mastitis events
newdata <- subset(uData, uData$health_event_mastitis_all != "NULL")
variable<-as.integer(paste(newdata$health_event_mastitis_all))
ggplot(newdata, aes(x = variable)) +
geom_bar(colour= "black", fill = mycoloursP[12]) + xlab("Mastitis events")+
theme_classic(20)+
scale_x_continuous(breaks=c(0:8))
#lameness events
newdata <- subset(uData, uData$health_event_lame_all != "NULL")
variable<-as.integer(paste(newdata$health_event_lame_all))
ggplot(newdata, aes(x = variable)) +
geom_bar(colour= "black", fill = mycoloursP[13]) + xlab("Lameness events")+
theme_classic(20)+
scale_x_continuous(breaks=c(0:11))
saving data
sessionInfo()## R version 3.6.2 (2019-12-12)
## Platform: x86_64-apple-darwin15.6.0 (64-bit)
## Running under: macOS Catalina 10.15.5
##
## Matrix products: default
## BLAS: /Library/Frameworks/R.framework/Versions/3.6/Resources/lib/libRblas.0.dylib
## LAPACK: /Library/Frameworks/R.framework/Versions/3.6/Resources/lib/libRlapack.dylib
##
## locale:
## [1] en_GB.UTF-8/en_GB.UTF-8/en_GB.UTF-8/C/en_GB.UTF-8/en_GB.UTF-8
##
## attached base packages:
## [1] grid stats graphics grDevices utils datasets methods
## [8] base
##
## other attached packages:
## [1] ggsci_2.9 survminer_0.4.7 ggpubr_0.4.0 survival_3.2-3
## [5] dplyr_1.0.0 plyr_1.8.6 ggthemes_4.2.0 magrittr_1.5
## [9] stringr_1.4.0 data.table_1.12.8 mgcv_1.8-31 gridExtra_2.3
## [13] lmerTest_3.1-2 lme4_1.1-23 Matrix_1.2-18 nlme_3.1-148
## [17] ggplot2_3.3.2
##
## loaded via a namespace (and not attached):
## [1] Rcpp_1.0.5 lattice_0.20-41 tidyr_1.1.0
## [4] zoo_1.8-8 digest_0.6.25 R6_2.4.1
## [7] cellranger_1.1.0 backports_1.1.8 evaluate_0.14
## [10] highr_0.8 pillar_1.4.6 rlang_0.4.7
## [13] curl_4.3 readxl_1.3.1 minqa_1.2.4
## [16] car_3.0-8 nloptr_1.2.2.2 rmarkdown_2.3
## [19] labeling_0.3 splines_3.6.2 statmod_1.4.34
## [22] foreign_0.8-76 munsell_0.5.0 broom_0.7.0
## [25] compiler_3.6.2 numDeriv_2016.8-1.1 xfun_0.15
## [28] pkgconfig_2.0.3 htmltools_0.5.0 tidyselect_1.1.0
## [31] tibble_3.0.3 km.ci_0.5-2 rio_0.5.16
## [34] crayon_1.3.4 withr_2.2.0 MASS_7.3-51.6
## [37] xtable_1.8-4 gtable_0.3.0 lifecycle_0.2.0
## [40] KMsurv_0.1-5 scales_1.1.1 zip_2.0.4
## [43] stringi_1.4.6 carData_3.0-4 farver_2.0.3
## [46] ggsignif_0.6.0 ellipsis_0.3.1 survMisc_0.5.5
## [49] generics_0.0.2 vctrs_0.3.1 boot_1.3-25
## [52] openxlsx_4.1.5 tools_3.6.2 forcats_0.5.0
## [55] glue_1.4.1 purrr_0.3.4 hms_0.5.3
## [58] abind_1.4-5 yaml_2.2.1 colorspace_1.4-1
## [61] rstatix_0.6.0 knitr_1.29 haven_2.3.1